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PREFACE 


The  project  on  which  we  report  here  had  three  objectives.  The  first  was  to  develop  new 
methods  for  modeling  regional  wave  propagation  in  laterally  heterogeneous  media.  Toward 
this  end  we  developed  some  new  computational  techniques  for  seismic  wave  scattering  based 
on  the  boundary  element  method  and  multipole  expansions.  In  this  report  we  include  a 
preprint  of  a  paper  entitled,  “Multiple  multipole  expansions  for  elastic  scattering,”  which 
has  been  submitted  to  Journal  of  the  Acoustical  Society  of  America.  The  second  objective 
was  to  develop  techniques  for  relative  event  location  and  apply  them  to  problems  of  nuclear 
monitoring.  The  second  part  of  this  report  is  the  application  of  a  relative  event  location 
method  we  developed,  based  on  waveform  cross-correlation,  to  the  location  of  presumed  ex¬ 
plosions  at  the  Balapan,  Kazakhstan  Test  Site.  Events  at  this  test  site  have  been  located 
very  accurately  on  the  basis  of  satellite  images,  and  this  provides  the  opportunity  to  test 
the  accuracy  of  our  approach  and  compare  it  to  the  location  accuracy  achievable  with  con¬ 
ventional  seismic  techniques.  Our  third  objective  was  to  study  the  effects  of  non-spherical 
cavities  on  the  seismic  wave  radiation  from  explosions.  The  final  section  of  this  report  is  a 
preprint  of  a  paper,  “Multipole  radiation  of  seismic  waves  from  explosions  in  non-spherical 
cavities,”  which  will  appear  in  the  Journal  of  the  Acoustical  Society  of  America. 
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MULTIPLE  MULTIPOLE  EXPANSIONS  FOR  ELASTIC 

SCATTERING 


Matthias  G.  Imhof 
Earth  Resources  Laboratory 

Department  of  Earth,  Atmospheric,  and  Planetary  Sciences 
Massachusetts  Institute  of  Technology 
42  Carleton  Street 
Cambridge,  MA  02142  -  1324 


ABSTRACT 

The  paper  presents  a  new  approach  to  solve  scattering  of  elastic  waves  in  two  dimensions.  Tradi¬ 
tionally,  wave  fields  are  expanded  into  an  orthogonal  set  of  basis  functions.  Unfortunately,  these 
expansions  converge  rather  slowly  for  complex  geometries.  The  new  approach  enhances  convergence 
by  summing  multiple  expansions  with  different  centers  of  expansions.  This  allows  irregularities  of 
the  boundary  to  be  resolved  locally  from  the  neighboring  center  of  expansion.  Mathematically,  the 
wavefields  are  expanded  into  a  set  of  non-orthogonal  basis  functions.  The  incident  wavefield  and 
the  fields  induced  by  the  scatterers  are  matched  by  evaluating  the  boundary  conditions  at  discrete 
matching  points  along  the  domain  boundaries.  Due  to  the  non-orthogonal  expansions,  more  match¬ 
ing  points  are  used  than  actually  needed  resulting  in  an  overdetermined  system  which  is  solved  in 
the  least  squares  sense. 

Since  there  are  free  parameters  such  as  location  and  number  of  expansion  centers  as  well  as 
kind  and  orders  of  expansion  functions  used,  numerical  experiments  are  performed  to  measure  the 
performance  of  different  discretizations.  An  empirical  set  of  rules  governing  the  choice  of  these 
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parameters  is  found  from  these  experiments.  The  resulting  algorithm  is  a  very  general  tool  to  solve 
relatively  large  and  complex  two-dimensional  scattering  problems. 

INTRODUCTION 

The  calculation  of  synthetic  seismograms  has  been  of  interest  for  many  years.  Various  methods 
have  been  proposed  for  modeling  waves  in  heterogeneous  media.  Each  of  them  has  it’s  own  range 
of  validity  and  interest.  Fully  numerical  techniques  in  the  space-time  domain,  either  in  the  finite 
difference  formulation  (Kelly  et  a/.,  1976;  Virieux,  1986)  (FD)  or  in  the  finite  element  formulation 
(Smith,  1974)  (FE),  handle  any  kind  of  waves  in  complex  media.  Unfortunately,  they  are  limited 
due  to  computer  memory  and  runtime  considerations,  since  for  many  problems  the  distances  be¬ 
tween  scatterer,  source  and  receiver  are  large.  An  area  containing  the  source,  the  receiver  and 
the  scatterers  plus  a  substantial  neighborhood  around  them  has  to  be  discretized,  which  might 
result  in  prohibitive  computation  times.  The  (generalized)  ray  theory  (Cerveny  et  a/.,  1977;  Hong 
and  Helmberger,  1978;  Cerveny  and  Psencik,  1984)  can  be  used  when  the  scatterers  and  their 
radii  of  curvature  are  large  compared  to  the  wavelength.  For  small  or  weak  scatterers,  the  Born 
approximation  (Miles,  1960)  allows  an  efficient  calculation  of  the  seismogram. 

In  other  cases,  the  problem  can  be  simplified  by  assuming  the  medium  to  consist  of  homogeneous 
regions  with  sharp  boundaries  inbetween.  Then,  reflectivity  (Muller,  1985;  Kennett,  1983)  and 
global  matrix  methods  (Chin  et  al,  1984)  are  routinely  used  for  planarly  or  cylindrically  layered 
media.  For  laterally  heterogeneous  media,  numerical  integration  over  wavenumber  can  be  used 
(Aki  and  Lamer,  1970;  Bouchon  and  Aid,  1977;  Haartsen  et  a/.,  1994).  In  the  elastic  case,  the 
classical  eigenfunction  expansion  (Morse  and  Feshbach,  1953)  (SMP)  allows  the  analysis  only  of 
simple  shapes  such  as  circular  cylinders  or  spheres  since  the  eigenmodes  do  not  decouple  otherwise 
(Pao  and  Mow,  1973).  Methods  based  on  the  perturbation  of  a  prescribed  geometry,  such  as  the 
T-matrix  method  (Waterman,  1976;  Bostrom,  1980)  work  extremely  well  for  certain  geometries  but 
are  harder  to  apply  efficiently  in  general  situations. 

The  method  we  present  is  a  derivative  of  the  boundary  element  methods  (Brebbia  and  Dominguez, 
1989)  (BEM).  It  was  first  presented  as  a  more  general  approach  for  electromagnetic  scattering 
(Ballisti  and  Hafner,  1983;  Hafner,  1990)  and  later  adapted  to  acoustic  scattering  problems  (Imhof, 
1995).  In  contrast  to  more  traditional  approaches,  the  wavefields  are  expanded  into  a  set  of  non- 
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orthogonal  and  non-complete  basis-functions.  Actually,  non-complete  basis-functions  are  not  a 
new  concept  since  for  numerical  and  computational  reasons  we  can  never  use  infinitely  many  basis- 
functions.  But  the  application  of  a  non-orthogonal  expansion  allows  the  reduction  of  the  truncation 
errors  (Hafner,  1993). 

To  solve  for  the  unknown  weighting  coefficients  of  the  basis-functions,  discrete  matching  points 
are  chosen  along  the  boundary  of  the  scattering  object.  In  the  elastic  case,  each  matching  point 
provides  four  boundary  conditions  and  thus  four  equations  involving  the  unknowns.  Because  the 
expansion  is  non-orthogonal,  we  require  more  equations  than  unknowns,  build  an  overdetermined 
matrix  system  and  solve  it  the  least-squares  sense.  Mathematically  speaking,  we  search  for  the  set 
of  weighting  coefficients  which  solves  the  problem  at  hand  “best”  employing  the  expansions  chosen. 
There  will  always  be  an  error  in  the  boundary  conditions  at  each  matching  point,  though  on  the 
average  these  errors  are  small.  Furthermore,  the  fields  in  between  matching  points  are  forced  to  be 
smooth,  such  that  no  wild  jumps  or  oscillations  can  occur.  Thus,  as  an  added  bonus,  we  control  the 
behavior  of  the  expansions  in  between  matching  points  where  we  have  no  control  using  traditional 
methods.  Also,  this  allows  us  to  see  in  which  parts  of  the  boundary  the  expansions  chosen  can 
solve  the  problem  and  where  they  need  further  refinement. 

This  paper  will  be  structured  as  follows:  First,  we  will  adapt  the  method  from  acoustical  (Imhof, 
1995)  to  elastic  scattering.  Then  we  present  results  from  several  calculations  and  compare  them  to 
solutions  obtained  by  the  finite  difference  method  and  the  classical  eigenfunction  expansion.  We 
show  how  different  discretizations  affect  the  solution.  Finally,  we  will  compile  these  findings  into 
an  empirical  set  of  rules  which  allows  us  to  set  a  problem  up  in  a  fashion  which  yields  satisfactory 
results  without  having  to  do  it  on  a  trial  and  error  basis. 

THEORETICAL  BACKGROUND 

We  would  like  to  model  how  an  incident  wavefield  u  *”‘^(x,a;)  scatters  from  an  object.  The  situation 
is  depicted  in  Figure  1.  For  the  sake  of  clarity,  we  will  suppress  the  time  factor  e”*"*  in  all  following 
expressions.  Superscripts  denote  the  region  to  which  a  material  property  or  fields  belongs  to,  and, 
to  distinguish  different  regions  or  domains,  we  use  the  symbol  F'^.  The  boundary  between  the  two 
regions  F®  and  F^  will  be  denoted  by  5Foi.  We  also  call  the  region  F°  the  background  and  define 
F^  to  be  the  scatterer. 
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In  the  frequency  domain,  elastic  P-SV  waves  travelling  in  a  two-dimensional,  homogeneous 
medium  are  described  by  (Pao  and  Mow,  1973) 


^  VV  •  u  —  xV  X  u  +  u  =  0 


(1) 


where  we  defined  the  wave  vectors  k  =  uj/ a  and  I  =  lo/P  for  a  particular  frequency  u  and  the 
propagation  velocities  a  =  \/X  +  2]i/p  and  /?  —  vWp-  The  parameter  p,  A  and  denote  the 
density  and  the  Lame  parameters  of  the  medium,  respectively. 

In  a  local  cylindrical  coordinate  system  (r,  0,^/)  centered  at  a  point  Xp  (Figure  2),  the  strains 
due  to  a  displacement  u  are  expressed  as  (Pao  and  Mow,  1973) 

(2a) 
(2b) 

v) 

All  other  components  are  zero  since  they  involve  the  Uy  component  or  cross-derivatives  with  respect 
to  y.  Thus,  the  stresses  are  linearly  related  to  the  strains  by  (Pao  and  Mow,  1973) 


dUr 

^rr  — 

dr 

duo  dUr 

^ee  = 

1/1  dUr  du0 

€r0  “  €.0Tr  ^ 

2\rdd~'^  dr 

^pq  —  X8pq  ^  ^  ^kk 
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where  Pj  ^  ^ 


(3) 


A  displacement  field  u  *^^(x)  incident  on  the  scatterer  will  induce  two  scattered  fields:  u^{x,u) 
outside  the  scattering  object  and  u^{x,uj)  on  the  inside.  The  displacements  and  stresses  inside 
and  outside  the  scatterer  are  related  by  the  boundary  conditions.  For  the  problem  posed,  these 
conditions  are  continuities  of  displacement  and  stresses  in  both  normal  and  tangential  direction. 
We  define  the  normal  h  to  point  from  medium  P^  into  medium  F^  as  depicted  in  Figure  2.  Using 
the  subscripts  n  and  t  to  denote  the  normal  and  tangential  direction,  we  write 


<+«r=“l 

1 


0  _|_ 

Ti  n.  ‘  ^  n  n. 


cr: 


+  ^ 


^  Tit 


nt 


(4a) 

(4b) 

(4c) 

(4d) 


Since  we  express  the  displacements  and  stresses  in  a  local  cylindrical  coordinate  system  {r,0,y), 
but  want  to  specify  the  boundary  in  a  local  cartesian  system  {n,t,y),  we  have  use  the  rotation 
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matrix  M  to  transform  the  individual  components 


U 


.  nty 


=  M  •  <T  •  M 


where  the  rotation  matrix  M  is  defined  by  the  unit  vectors  f,  6,  n  and  i. 

(  f).-r  n  - 6  0  \ 


M  = 


n-T  n.' 

i  ■  f  i  ■  9  0 
0 


0  1  ) 

Instead  of  using  the  displacement  u  directly,  we  break  it  into  two  parts 


(5) 

(6) 


(7) 


u{x)  ~  V$(f)  +  V  X  (8) 

using  the  scalar  potentials  $(x)  and  Then,  equation  (1)  separates  into  two  independent 

Helmholtz  equations: 

(V2  +  jfc2)$(£^^)  =  0  (9a) 

(V2+/2)'^-(f,a;)  =0  (9b) 

Therefore,  we  replace  the  induced  displacement  fields  u^{x^lo)  and  u^{x,uj)  by  the  potentials 
$°(x, O'),  ^'®(x, O’),  $^(x,  w)  and  w).  Similar  to  the  acoustic  case  (Imhof,  1995),  we  expand 

the  potential  fields  as: 


P  +A7 


=  E  E  +4> 

(10a) 

p=:0  n=:~N 

P 

(10b) 

p=0n=-Ar 


where  (l>pn{x^  x^^  u)  and  x^^  cj)  are  solutions  to  either  Helmholtz  equation  (9a)  or  (9b), 

respectively.  The  error  terms  and  are  included  not  only  since  the  series  are  truncated  after 
±N  terms  but  also  because  an  expansion  of  this  form  is  mathematically  non- orthogonal 

An  expansion  of  the  form  (10a)  or  (10b)  is  known  as  a  multiple  multipole  (MMP)  expansion 
(Ballisti  and  Hafner,  1983;  Hafner,  1990).  Setting  P  to  zero  yields  the  classical  eigenfunction 
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(SMP)  expansion  (Morse  and  Feshbach,  1953).  In  the  background  region  F^,  we  choose  propagatory 
solutions  for  the  expansion  functions  and 

=  H^n\  {k^  -  ^S\)  (11a) 

fppni^,  w)  =  {1°  \x  -  f  “D  e"”®  (lib) 

For  the  sake  of  clarity,  the  notation  for  an  expansion  center  means  the  center  for  the 
expansions  of  u)  and  The  location  is  not  fixed  yet.  We  can  either  set  the  center  Xp 

inside  the  scatterer  F^  and  thus  Xp  G  F^  or  we  can  place  it  in  the  background  {Xp  G  F^).  For  a  finite 
scatterer  F^,  we  have  two  possible  choices  for  the  expansions  of  the  fields  cj)  and  '^^{x^u). 
First,  we  can  place  the  expansion  centers  into  the  scatterer  itself  {Xp  G  F^)  and  use  expansions 
involving  the  Bessel  solutions  corresponding  to  standing  waves  (Morse  and  Feshbach,  1953). 


^pn(^>^p^^^^^)  =  {k^  l^-^pl)  if  ^p  e  r’-  (12a) 

^pn(^^^p^^^w)  =  (^M^-^pl)  iffp€r^  (12b) 

Second,  we  can  place  the  the  expansion  centers  into  the  background  {Xp  G  F^)  and  use  propagatory 
solutions  involving  the  Hankel  functions  of  the  first  kind. 

ln{x, X^,k^,uj)  =  lx  -  Xp^ I)  e*”®  if  x^  €  T®  (13a) 

V’p„(5, Xp,k\oj)  =  (/^  |x  -  Xp^ I)  if  Xp^  G  T®  (13b) 


These  expansions  represent  waves  propagating  from  the  expansion  center  toward  the  scatterer 
(Morse  and  Feshbach,  1953).  In  the  scatterer,  we  need  waves  propagating  in  all  the  directions. 
Thus,  we  place  expansion  centers  all  around  the  scatterer  and  “illuminate”  the  region  F^  from  all 
sides. 

To  emphasize  the  difference  between  the  expansions  (12)  and  (13);  if  Xp  is  placed  inside  the 
scatterer  F^,  then  we  have  to  use  (12)  because  the  Bessel  solutions  J|„|  e*”®  represent  standing 
waves.  Contrary  to  the  Hankel  solutions,  the  Bessel  solutions  do  not  have  a  singularity  at  their 
origin.  Expansion  functions  involving  Hankel  functions  may  never  be  used  for  the  wavefields 
of  the  domain  in  which  their  expansion  center  Xp  is  located  in  since  the  singularities  at  the  origin 
represent  sources.  But  by  definition,  the  only  source  in  the  problem  posed  is  the  incident  field 
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But  if  Xp  is  located  outside  the  scatterer,  then  we  have  to  use  (13)  with  the  Hankel  solutions 
^|n|^  representing  wavefields  propagating  towards  the  scatterer.  The  singularities  pose  no 
problem  anymore  since  they  are  not  located  in  the  domain  Figure  3  illustrates  this  subtlety. 

We  solve  for  the  unknown  coefficients  and  by  enforcing  the  boundary  conditions  (4a) 
-  (4d)  on  M  discrete  matching  points  Xm  along  the  domain  boundary  ^Foi.  Since  we  have  four 
boundary  conditions,  each  matching  point  also  provides  four  rows  of  the  linear  matrix  system. 
Altogether,  we  have  4 J  =  2  •  2  •  P  •  (2^  +  1)  unknown  coefficients  6^^.  To  simplify  the  no¬ 
tation,  we  eliminate  an  index  by  sequentially  renumbering  the  double-indexed  expansion  functions 
(f)p^{x^Xp,k^^u)  and  the  coefficients  which  results  in  (j)^{x^xf^k^^u)  and  a^,  respectively.  The 
same  is  done  with  the  'ipp^{x^Xp,l^,u)  and  the  coefficients  resulting  in  il)j{x^xf,l^,ijo)  and 
respectively.  Putting  all  together,  we  have  to  solve  a  matrix  system  of  the  form 


( 
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(14) 


4M 


where  we  used  the  submatrices  and  to  denote  the  normal  displacements  Un  at  the  matching 
points  due  to  and  respectively.  The  submatrices  and  are  the  same  but  for  the 
tangential  displacements  The  submatrices  and  contain  the  normal  stresses  anm  while 
and  contain  the  tangential  stresses  ant- 

Defining  the  matching  points  by  their  location  Xm^  we  can  write  these  submatrices  as 


II 

^n,mj  ~ 

(15a) 

^t,mj  ^  i^m)) 

(15b) 

^nn,mj  ~  ^nn(<^j  (^m)) 

^nn,Tnj  ~  ^nni^P j{^m)) 

(15c) 

b 

It 

(15d) 

where  we  used  the  index  m  G  ,M}  to  denote  the  matching  points  the  index  j  G 

,  J}  for  the  expansion  functions  (j)j^  'ipj  and  the  index  d  G  {0, 1}  for  the  domain.  The  notation 
Un{(pj{xm))  stands  for  the  normal  displacement  due  to  the  expansion  function  (j)j  evaluated  at  the 
matching  point  Xm-  The  other  ones  are  to  be  interpreted  similarly. 
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The  vectors  a^,  6^,  and  in  equation  (14)  contain  the  unknown  coefficients  a^,  6],  aj,  for 
the  expansion  functions  '0^,  (f)^  and  respectively.  The  vectors  Ut,  ann,  and  ant  hold  the 
normal  and  tangential  displacements  as  well  as  normal  and  tangential  stress  at  the  M  matching 
points  due  to  the  incident  field  t? 

ut,m  =  (16a) 

=  <T(^m)  (16b) 

Finally,  the  matrix  equation  (14)  contains  the  residual  vectors  e^,  e^,  e^nj  ^nt  with  the  misfit  of  the 
boundary  conditions  at  the  individual  matching  points. 

NUMERICAL  RESULTS 

To  reduce  numerical  noise,  we  make  the  materials  slightly  lossy  by  adding  a  small  imaginary 
component  ui  to  the  frequency  (Bouchon  and  Aki,  1977).  Thus,  we  have  to  evaluate  Bessel  functions 
with  a  complex  argument  (Amos,  1986).  After  the  transformation  from  the  frequency  domain  into 
the  time  domain,  we  recover  the  true  amplitude  by  a  multiplication  with  The  matrix  system 
is  solved  by  QR  decomposition  using  Givens  rotations  (Wilkinson,  1988)  which  allows  to  build 
the  matrix  system  row  by  row  while  only  a  triangular  matrix  with  dimensions  of  the  number  of 
unknowns  has  to  be  kept  in  memory  (George  and  Heath,  1980).  Since  we  want  to  calculate  synthetic 
seismograms  using  a  frequency  domain  method,  we  have  to  solve  the  scattering  problem  for  a  range 
of  frequencies  and  later  apply  a  Fourier  transformation  to  obtain  the  seismograms.  All  these 
problems  can  be  solved  independently  of  each  other.  Consequently,  the  algorithm  is  implemented 
on  an  nCUBE2  parallel  computer  where  each  processor  will  calculate  a  few  frequencies. 

We  will  now  show  how  the  method  performs  solving  a  very  simple  problem  using  different 
ways  to  discretize  it.  For  the  sake  of  simplicity,  the  incident  field  u  is  an  explosive  line  source 
modulated  with  a  Ricker  wavelet  (Hosken,  1988;  Paillet  and  Cheng,  1991)  of  50  Hz  center  frequency. 
Altogether,  64  receivers  will  measure  the  Uz  component  of  the  scattered  field  it  ^  in  the  background. 

The  rather  generic  scatterer  is  depicted  in  Figure  4.  Its  size  is  roughly  240m  in  length  and 
50m  thickness.  The  velocities  in  the  background  are  =  2000m/s  and  =  1155m/s,  while  in 
the  scatterer  they  are  =  3000m/s  and  =  1732m/s.  Thus,  the  Poisson’s  ratio  is  the  same  for 
both  regions  (cr  =  0.25).  To  facilitate  the  comparison  with  a  solution  obtained  by  finite  differences. 


Un,m  = 

ann^m  —  ^nn  \^m) 
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the  density  p  =  =  pMs  kept  constant  at  2000kg/m^.  The  center  frequency  of  50Hz  yields 

wavelengths  of  essentially  of  same  size  as  the  scatterer. 

In  order  to  have  a  reference  seismogram  to  compare  the  different  solutions  with,  we  calculate 
the  solution  using  a  finite  difference  (FD)  method  (Kelly  et  al.^  1976;  Peng  and  Toksoz,  1994). 
The  resulting  reference  seismogram  is  shown  in  Figure  5.  As  a  measure  of  how  well  the  MMP 
seismogram  t)  correlates  with  the  FD  reference  seismogram  t),  we  define  the  root 

mean  square  error  RMS  by  summing  over  the  squared  difference  between  the  two  seismograms 


RMS  = 


R  T 


EE{  uMMP(^r,t)  -  uf^{r,t)y 


(17) 


where  Ur{T^t)  denotes  the  vertical  displacement  measured  at  recorder  r  at  time  sample  t.  it!  =  64 
is  the  number  of  recorders  and  T  —  256  is  the  total  number  of  time  samples. 


MMP  versus  the  Finite  Difference  Reference  Solution 

As  a  first  example,  we  show  both  a  solution  obtained  by  MMP  expansions  and  the  reference  solution 
as  obtained  by  finite  differences.  For  the  finite  difference  case,  we  used  a  grid  spacing  of  Im  and 
a  grid  of  750  by  750  points.  The  grid  dimensions  are  larger  than  needed  to  prevent  any  reflection 
from  the  boundaries  to  reach  the  receivers.  The  timestep  was  0.05ms.  The  runtime  on  a  nCUBE2 
using  64  nodes  was  23  minutes.  The  seismogram  calculated  by  finite  differences  is  presented  in 
Figure  5. 

For  the  MMP  expansion,  we  used  a  total  of  256  expansion  functions,  128  matching  points,  8 
expansion  centers  and  64  frequencies.  The  resulting  runtime  on  a  nCUBE2  using  again  64  nodes 
was  12  minutes.  The  two  methods  yield  very  similar  results.  Figure  6  shows  the  seismogram 
calculated  using  the  MMP  expansion.  As  can  be  seen,  they  agree  very  well  in  both  traveltimes  and 
phases. 


Effect  of  the  Number  of  Expansion  Functions 

As  a  second  experiment,  we  study  how  the  number  of  expansion  functions  affects  the  solutions 
obtained.  We  start  with  totally  32  expansion  functions  located  at  8  expansion  centers.  Thus,  we 
have  one  monopole  for  each  potential  and  each  region  at  every  expansion  center. 
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We  calculate  the  resulting  seismogram  and  estimate  the  RMS  error.  Then  we  double  the  number 
of  expansion  functions  per  expansion  center,  calculate  the  seismograms  anew,  estimate  the  RMS 
error  and  so  on  until  1024  expansion  functions  are  used.  The  number  of  matching  points  is  kept 
constant  at  M  =  1024  while  the  number  of  expansion  centers  is  kept  constant  at  P  =  8.  Figure  7 
shows  the  resulting  RMS  error  as  a  function  of  the  total  number  of  expansion  functions  used.  A 
first  observation  is  that  256  expansion  functions  seems  to  be  the  critical  amount.  Using  fewer  yields 
solutions  which  cannot  capture  important  features  of  the  reference  seismogram,  the  solutions  do 
not  converge.  Figure  8  shows  a  seismogram  which  is  typical  for  a  not  converged  solution.  The 
seismogram  was  obtained  with  64  expansion  functions.  For  more  than  256  expansion  functions,  we 
have  convergence  where  RMS  error  decreases  slowly  with  increasing  number  of  expansions  functions. 


MMP  versus  SMP  Expansion 


The  next  numerical  experiment  we  perform  is  to  show  the  enhanced  convergence  of  the  MMP 
expansion  compared  to  the  classical  eigenfunction  expansion  (SMP).  As  mentioned  priorly,  the 
eigenfunction  expansion  corresponds  to  an  expansion  (10a)  or  (10b)  with  only  one  expansion  center. 
Thus,  we  perform  the  same  experiment  as  before  but  use  only  one  expansion  center.  Again,  we  start 
with  one  expansion  function  per  domain  and  potential  which  yields  totally  4  expansion  functions. 
We  calculate  the  resulting  seismograms  and  estimate  the  RMS  error.  The  resulting  seismogram 
is  presented  in  Figure  9.  The  seismogram  is  clean  enough  to  be  mistaken  as  correct  but  has  no 
resemblance  with  the  reference  solution  shown  in  Figure  5.  Then  we  double  the  number  of  expansion 
functions  per  expansion  center,  calculate  the  seismograms,  estimate  the  RMS  error  and  so  on.  The 
number  of  matching  points  is  kept  constant  at  M  =  1024. 

Figure  7  shows  the  resulting  RMS  error  as  a  function  of  the  total  number  of  expansion  functions 
used.  We  notice  that  the  MMP  expansion  using  8  expansion  centers  always  performs  better. 
Unfortunately,  using  more  than  256  expansion  functions  in  the  SMP  expansion  yields  no  useful 
result  anymore  because  the  expansion  functions  of  higher  order  violate  the  sampling  condition 
(Hafner,  1990;  Imhof,  1995).  The  maximum  order  of  a  multipole  is  given  by  the  largest  angle 
^max  between  any  two  adjacent  matching  points  and  the  location  of  the  multipole: 

TT 


< 


(18) 


The  increased  error  in  the  SMP  expansion  for  32  and  64  expansion  function  is  an  effect  of  the 
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error  measure  (17)  which  cannot  acount  for  phase  shifts.  Contrary  to  the  MMP  expansion,  a  SMP 
expansion  cannot  solve  the  problem  posed  in  Figure  5. 

Effect  of  Number  and  Location  of  Expansion  Centers 

The  next  numerical  experiment  is  to  examine  the  importance  and  effect  of  the  number,  location  and 
distribution  of  expansion  centers.  As  priorly  mentioned,  we  have  the  choice  of  placing  the  expansion 
centers  for  the  and  fields  either  in-  or  outside  the  scatterer  and  thus  expanding 

either  into  standing  waves  J\n\  {kr)  or  into  propagating  waves  (kr)  respectively.  We 
will  use  both  to  study  the  difference. 

We  calculate  the  solutions  for  a  range  of  expansion  centers  while  keeping  the  total  number  of 
expansion  functions  constant  at  256.  Also,  the  number  of  matching  points  is  kept  constant  at  256. 
The  overall  computational  effort  to  calculate  one  seismogram  is  kept  constant.  The  resulting  RMS 
errors  are  shown  in  Figure  10.  It  is  surprising  how  broad  the  ‘U’  shaped  minimal-error  region  is. 
The  range  from  4  up  to  17  expansion  centers  converges.  Indeed,  the  minimal  RMS  error  obtained 
by  11  expansion  center  is  only  slightly  better  any  other  discretization  employing  4  to  17  centers. 
Remarkably,  MMP  expansions  seem  to  be  very  insensitive  to  the  discretization  used!  Neither  the 
number  of  expansion  centers  nor  the  kind  of  expansions  changes  the  RMS  errors  by  much,  although 
the  use  of  (kr)  produces  a  smoother  RMS  error  curve. 

The  pathological  case  with  23  expansion  centers  shows  that  the  RMS  error  finally  increases 
when  more  and  more  expansion  centers  are  used.  In  this  particular  case,  the  expansion  centers 
were  separated  by  only  a  quarter  of  the  dominant  wavelength.  The  different  expansion  functions 
begin  to  interact  by  approximating  higher  order  solutions  to  the  wave  equation:  it  is  well  known 
that  two  monopoles  of  opposite  sign  placed  closely  together  are  equivalent  to  a  dipole.  Thus,  the 
matrix  system  becomes  more  and  more  ill-conditioned  since  each  expansion  center  could  be  replaced 
by  the  adjacent  ones.  Moreover,  we  add  more  and  more  similar  equations  to  the  matrix  system 
which  renders  it  more  and  more  ill-conditioned. 

For  comparison,  we  also  use  a  simple  boundary  element  (BEM)  discretization  with  the  same 
number  of  matching  points  and  expansion  functions.  Along  the  boundary  inbetween  matching 
points,  we  place  rotational  and  compressional  monopole  sources.  As  in  the  MMP  cases,  we  use 
point  matching  and  solve  the  system  in  the  least-squares  sense.  The  resulting  large  RMS  error 
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indicates  that  the  seismogram  obtained  is  not  correct.  Indeed,  it  contains  mainly  the  reflections 
from  the  top  of  the  scatterer.  Both  reflections  from  the  bottom  and  internal  multiple  scattering 
are  mostly  missing. 

Effect  of  the  Number  of  Matching  Points 

The  last  numerical  experiment  examines  how  the  number  of  matching  points  affects  the  solutions. 
Actually,  not  the  number  of  matching  points  but  the  ratio  between  the  total  number  of  equations 
in  the  matrix  system  and  the  number  of  expansion  functions  used  is  the  important  parameter. 
In  accordance  with  the  earlier  experiments,  we  choose  12  expansion  centers  and  keep  the  number 
of  expansion  functions  constant  at  256.  Since  each  matching  point  provides  4  equations  (one  for 
each  boundary  condition),  we  start  out  with  64  matching  points  along  the  boundary  which  provide 
256  equations  altogether.  We  calculate  the  resulting  seismogram,  estimate  the  RMS  error,  double 
the  number  of  matching  points  and  so  on.  Figure  11  shows  the  RMS  error  against  the  number 
of  equations  per  expansion  function.  Since  the  expansion  is  non-orthogonal,  it  is  not  surprising 
that  we  get  a  large  RMS  error  when  we  use  as  many  equations  as  we  have  unknowns.  Using  twice 
as  many  equations  as  unknowns  provides  the  optimal  result.  Afterwards,  the  more  equations  we 
add,  the  more  the  RMS  increases  since  the  matrix  system  becomes  more  ill-conditioned  with  each 
additional  equation  we  add.  The  result  is  more  errors  due  to  roundoff  and  other  numerical  effects. 

Using  twice  as  many  equations  as  unknowns  yields  a  distance  of  4m  between  matching  points. 
This  spacing  corresponds  to  10  matching  points  per  dominant  wavelenght  (40m).  Assuming  that 
the  highest  frequency  in  the  propagating  seismic  wavelet  is  3  times  the  center  frequency  of  50Hz, 
the  boundary  is  sampled  with  3  matching  points  per  wavelength  for  the  highest  frequency.  The 
sampling  theorem  which  states  that  the  boundary  has  to  be  sampled  at  least  twice  per  wavelength  to 
prevent  aliasing  (Bouchon  and  Aki,  1977),  is  just  satisfied.  Thus,  it  is  also  theoretically  reasonable 
to  have  about  10  matching  points  per  dominant  wavelength. 

DISCUSSION 

Combining  these  numerical  experiments  with  prior  experiences  with  electromagnetic  (Hafner,  1990) 
and  acoustical  MMP  methods  (Imhof,  1995),  we  obtain  a  set  of  empirical  rules  how  to  discretize 
elastic  scattering  problems.  A  very  important  parameter  is  the  radius  of  greatest  influence  of  a 
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multipole  which  is  \/2  times  the  distance  between  the  center  of  expansion  and  the  closest  matching 
point, 

•  The  radius  of  greatest  influence  should  be  on  the  order  of  the  dominant  wavelength. 

•  No  expansion  center  should  be  within  the  radius  of  greatest  influence  of  any  other  expansion 
center. 

•  There  should  be  10  matching  points  per  dominant  wavelength 

•  There  has  to  be  at  least  half  a  matching  point  per  expansion  function  or  similarly  two  equa¬ 
tions  per  expansion  function. 

•  The  maximum  order  N  of  a.  multipole  is  given  by  the  sampling  theorem:  N  < 

•  Expansions  of  the  form  {kr)  should  not  be  used  for  the  region  their  expansion  center 
Xp  is  located  in. 

All  of  these  rules,  except  the  last  one,  are  only  general  guidelines.  Adhering  to  these  guidelines 
yields  satisfactory  results.  As  the  numerical  experiments  show,  all  parameters  can  be  varied  by 
large  amounts  while  only  perturbing  the  resulting  solution.  The  MMP  method  is  not  very  sensitive 
to  the  actual  discretization  used. 

As  shown,  the  MMP  expansion  converges  faster  than  the  classical  multipole  or  simple  boundary 
element  expansions  for  complex  scattering  geometries.  The  method  is  either  able  to  solve  scattering 
problems  involving  harmonic  sources  or  to  calculate  seismograms  by  Fourier  synthesis.  For  the 
problem  posed,  we  also  found  the  MMP  expansions  to  be  faster  than  finite  difference  modeling 
with  approximately  the  same  degree  of  accuracy.  In  the  example,  source,  receivers  and  scatterers 
were  located  close  to  each  other.  For  problems  with  larger  distances  between  them,  we  can  expect 
an  even  greater  decrease  in  computation  time  compared  to  finite  differences.  Furthermore,  due 
to  its  spectral  nature,  attenuation  can  easily  be  accounted  for.  Thus,  the  MMP  method  is  well 
suited  for  a  large  range  of  scattering  problems  since  both  acoustic  and  elastic  media  with  different 
boundary  conditions  (fluid-fluid,  elastic-elastic  and  others)  can  be  treated  exactly  the  same  way  in 
this  algorithm. 
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Figure  1:  Schematic  representation  of  the  scattering  experiment.  An  incident  field 
illuminates  a  bounded  two-dimensional  inhomogeneity  which  induces  a  scattered  field  a;)  in 

the  background  medium  as  well  as  a  field  u^{x,u)  in  the  scatterer  itself 


Figure  2:  Schematic  of  the  coordinate  systems  used.  Additional  to  a  global  cartesian  coordinate 
frame  local  cylindrical  systems  (r^O^y)  with  origins  at  Xp  are  used.  Such  a  local  origin  or 

expansion  center  is  depicted  by  the  triangle.  Boundaries  between  different  media  are  defined  by 
discrete  matching  points  located  at  Xm  where  normal  n  and  tangential  t  directions  are  specified. 
The  matching  points  are  denoted  by  squares. 
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Figure  3:  Basis  functions  can  either  contain  Hankel  functions  Hn  or  Bessel  functions  J„.  If  the 
same  expansion  center  is  to  be  used  for  as  for  (j)^,  then  the  inside  field  has  to  be  expanded 
using  the  Bessel  functions  Jn  since  they  represent  standing  waves.  If  the  inside  and  the  outside 
scattered  field  are  to  represented  by  Hankel  functions  Hn^  different  expansion  centers  have  to  be 
used.  Expansion  centers  are  depicted  by  a  triangle  while  squares  are  used  for  matching  points. 
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Figure  4:  Generic  scatterer  used  for  numerical  experiments.  The  scatterer  is  illuminated  by  an 
explosive  line  source  modulated  by  a  Ricker  wavelet  of  50Hz  center  frequency.  The  velocities  in  the 
background  domain  F*’  are  =  2000m/s  and  =  1155m/s.  The  velocities  in  the  scatterer  F^ 
are  =  3000m/s  and  =  1732m/s.  The  Poisson’s  ratio  is  the  same  for  both  regions  (a  =  0.25). 
For  simplicity,  the  density  is  kept  constant  at  p  =  2000kg/m®. 
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0.0 


Figure  6:  The  seismogram  of  the  model  shown  in  Figure  4  calculated  using  the  MMP  algorithm. 
Altogether,  256  expansion  functions,  8  expansion  centers,  128  matching  points  and  64  frequencies 
were  used.  As  can  be  seen,  the  MMP  solution  agrees  very  well  with  the  finite  difference  reference 
seismogram  shown  in  Figure  5. 
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0.3 


RMS  Error  vs  Total  Number  of  Expansion  Functions 


Figure  7:  Comparison  between  the  traditional  eigenfunction  expansion  SMP  (dashed)  and  the 
MMP  expansion  (solid).  Shown  is  how  the  total  number  of  expansion  functions  affects  the  RMS 
error  compared  to  the  FD  reference.  The  SMP  actually  never  converges  since  for  512  expansion 
functions  it  violates  the  sampling  condition.  For  256  and  more  expansion  functions,  the  MMP 
expansion  converge. 
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Figure  8:  The  seismogram  for  the  case  when  64  expansion  functions  are  used.  The  seismogram  is 
very  noisy.  Some  of  the  prominent  features  in  Figure  5  begin  to  show  up,  but  the  expansions  have 
not  converged  yet.  More  terms  have  to  be  used  to  obtain  convergence. 
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Figure  9:  The  seismogram  for  the  case  when  only  4  expansion  functions  are  used.  Clearly,  no  self¬ 
interaction  of  the  scattered  wavefields  is  possible.  Unfortunately,  the  seismogram  is  clean  enough 
to  be  mistaken  as  correct. 
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RMS  Error  vs.  Number  of  Expansion  Centers  (Total  Number  of  Expansions  =  256) 
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Figure  10:  The  effect  of  number  and  location  of  the  expansion  centers.  The  total  number  of 
expansion  functions  is  kept  constant  at  256  while  the  number  of  expansion  centers  is  varied  from  1 
up  to  23.  Expansions  using  the  same  expansion  centers  G  for  and  thus  Bessel 

functions  Jn  as  well  as  expansions  using  expansion  centers  G  for  'ijp  and  x^  G  F^  for  -0^ 
and  thus  Hankel  functions  Hn  are  tested.  The  difference  between  these  two  kinds  of  expansions  is 
rather  small.  Placing  all  expansion  centers  onto  the  boundary  and  using  only  the  0^^  order  terms 
which  corresponds  to  a  simple  boundary  element  expansions  fails  surprisingly. 
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RMS  Error  vs.  Number  of  Equations 


Figure  11:  Influence  of  the  number  of  matching  points  on  an  expansion  with  256  expansion  functions 
and  8  centers  of  expansion.  Each  matching  point  provides  4  equations.  Since  the  expansion  is  non- 
orthogonal,  using  as  many  equations  as  unknowns  to  be  resolved  does  not  yield  a  correct  result. 
Adding  more  and  more  equations  to  the  system  increases  the  condition  number  and  thus  the  RMS 
error  is  increased  due  to  numerical  errors. 
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SUMMARY 

We  apply  a  technique  for  relative  event  location  with  differential  arrival  times  to  relocate 
four  presumed  nuclear  explosions  occurring  at  the  Balapan,  Kazakhstan  test  site  during  1987. 
Waveform  cross-correlation  analysis  was  performed  on  seismograms  from  ten  far-regional  and 
teleseismic  stations,  to  determine  differential  arrival  times  for  P  waves  and  some  PcP  waves 
relative  to  the  Joint  Verification  Event  of  14  September  1988.  The  resulting  differential  arrival 
time  data  were  used  to  determine  the  epicenters  of  the  four  events  relative  to  the  epicenter  of 
the  JVE  event.  We  compare  our  locations  to  highly  precise  epicenters  determined  from  satellite 
images,  and  to  teleseismic  locations  obtained  with  the  conventional  location  method  applied  to 
hundreds  of  absolute  arrival  time  picks.  Our  locations  differ  from  the  satellite  locations  by  2  to 
3  kilometers  and  are  more  accurate  than  the  locations  obtained  from  hundreds  of  data. 

INTRODUCTION 

The  ability  to  accurately  locate  seismic  events  is  of  crucial  importance  for  seismic  identification 
and  discrimination.  We  have  developed  a  high  precision  relative  event  location  method  based 
on  seismic  waveform  cross-correlation  analysis.  Previously,  we  reported  the  application  of  this 
method  to  quarry  blasts  in  Estonia  using  near  regional  data  from  the  Scandinavian  arrays 
FINESA,  NORESS  and  ARCESS  (Toksoz  et  al^  1993;  Rodi  et  a/.,  1994).  The  results  of  thus 
study  showed  a  significant  improvement  in  location  accuracy  compared  to  routine  locations 
done  by  the  Intelligent  Monitoring  System,  with  relative  epicenters  between  events  estimated 
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to  have  accuracies  of  approximately  one  kilometer.  In  the  present  report,  we  extend  our  method 
to  locate  nuclear  explosions  in  the  Balapan,  Kazakhstan  test  site  (KTS)  using  far-regional  and 
teleseismic  waveform  data. 

During  the  years  1965  through  1989,  101  presumed  nuclear  explosions  from  the  Balapan, 
Kazakhstan  Test  Site  (KTS)  were  detected  and  located  teleseismically  (e.g,,  Marshall  et  a/., 
1984;  Ringdal  et  al,  1992).  There  have  been  a  number  of  teleseismic  location  studies  of  Balapan 
events  using  the  joint  epicenter  determination  (JED)  method  of  Douglas  (1967).  Marshall  et 
al  (1984)  located  61  Balapan  explosions  through  1982,  using  the  location  of  event  650115, 
which  is  determined  by  LANDSAT  satellite  imaging,  as  a  master  event.  Lilwall  and  Farthing 
(1990)  relocated  presumed  nuclear  explosions  from  1973  to  1989  using  seven  master  events. 
Furthermore,  Thurber  et  al.  (1994)  used  27  master  events  to  constrain  the  teleseismic  locations 
with  both  JED  and  a  master  event  location  algorithm.  The  locations  of  an  additional  20 
master  events  w^ere  determined  by  analyzing  a  SPOT  satellite  image  (Thurber  et  al..,  1993). 
The  locations  of  these  presumed  nuclear  explosions  have  also  been  routinely  determined  by  ISC 
and  NEIC.  The  ISC  and  NEIC  generally  used  several  hundred  arrival  time  data  recorded  by 
all  available  stations  around  the  world  to  determine  the  epicenters.  However,  none  of  these 
previous  location  studies  have  utilized  the  seismic  waveforms  recorded  from  the  events. 

Using  both  SPOT  and  LANDSAT  satellite  images,  Thurber  et  al.  (1993,  1994)  determined 
the  locations  of  101  presumed  explosions  at  KTS  with  100  m  to  200  m  precision.  These  studies 
found  that  the  teleseismic  locations  (JED)  of  most  events  agreed  with  the  LANDSAT/  SPOT 
shot  points  to  within  about  1  to  2  km.  The  high  precision  locations  based  on  satellite  images 
also  provide  a  ground  truth  database  for  testing  our  relative  location  method.  In  this  report, 
we  apply  our  method  to  seismic  waveform  data  from  five  nuclear  explosions  recorded  at  ten 
stations,  with  epicentral  distances  from  23  to  67  degrees. 

DATA  PROCESSING 

Figure  1  is  a  world  map  showing  locations  of  five  presumed  nuclear  explosions  (circles)  at 
Balapan,  KTS,  and  ten  seismic  stations  (triangles)  used  in  this  study.  All  stations  have  vertical 
short-period  seismic  instruments  and  a  few  of  them  also  were  equipped  with  three-component 
broadband  sensors.  The  sampling  rate  is  from  20  to  40  samples  per  second.  Figure  2  shows 
the  locations  of  the  five  presumed  explosions,  as  determined  by  satellite  imaging  (Thurber  et 
aZ.,  1993)  and  the  PDE  using  teleseismic  arrival  time  data  (NEIC/USGS,  1987,  1988).  The 
hypocentral  parameters  of  the  events  are  also  listed  in  Table  1.  In  order  to  examine  the  relative 
locations  among  the  events,  we  use  the  satellite  image  location  of  event  94,  the  JVE  (Joint 
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Table  1:  Hypocentral  Parameters  of  Five  Presumed  Underground  Nuclear  Explosions  at  the 
Balapan  Test  Site 


Event 

YrMoDa 

HrMinSec 

Latitude 

Longitude 

mb 

(a) 

(b) 

(b) 

(c) 

85 

870620 

005304.8 

49.9367 

78.7464 

6.0 

86 

870802 

005806.8 

49.8806 

78.8750 

5.8 

88 

871213 

032104.9 

49.9614 

78.7933 

6.1 

89 

871227 

030504.8 

49.8789 

78.7253 

6.0 

94 

880914 

035957.4 

49.8781 

78.8239 

6.0 

Notes:  (a)  Origin  time  is  from  NEIC/USGS.  (b)  Latitude  and  longitude  are  from  from  Thurber 
et  ai  (1993).  (c)  is  from  Ringdal  et  ai  (1992), 

Verification  Experiment)  event,  as  a  common  reference  point.  In  Figure  2,  the  locations  of  the 
other  four  events  are  plotted  relative  to  this  reference  point.  We  can  see  from  the  figure  that 
location  differences  between  the  satellite  image  and  PDE  locations  vary  from  2  to  10  km.  We 
will  use  our  waveform  correlation/relative  event  location  algorithm  to  examine  what  relative 
location  accuracy  can  be  achieved  with  waveforms  from  only  ten  stations.  The  stations  are 
listed  in  Table  2,  together  with  their  distances  and  azimuths  from  the  JVE  event. 

Figure  3  shows  the  seismograms  of  the  JVE  event  (mfo=6.0)  recorded  at  the  ten  stations 
shown  in  Figure  1.  We  took  the  JVE  event  to  be  a  master  event,  and  determined  differential 
arrival  times  between  it  and  the  other  four  events.  Figure  4  shows  vertical  seismograms  of 
the  five  study  events  recorded  at  station  LZH,  with  an  epicentral  distance  of  23  degrees  and 
an  azimuth  of  118  degrees.  We  note  that  the  P  waveforms  for  the  first  few  seconds  are  very 
similar  to  one  another,  implying  these  events  do  not  have  significant  hypocentral  separation. 
Since  the  epicentral  distance  from  station  LZH  to  the  test  site  is  at  a  far-regional  distance, 
the  P  waveforms  are  more  complex  than  those  recorded  at  teleseismic  distance.  This  can  be 
seen  in  Figure  5,  which  shows  vertical  component  seismograms  of  the  five  events  recorded  at 
station  COL,  at  an  epicentral  distance  of  60  degrees  and  an  azimuth  of  21  degrees.  In  contrast 
to  Figure  4,  the  P  waveforms  recorded  at  station  COL  are  relatively  simple.  Like  LZH,  the 
waveforms  of  the  five  events  are  similar  to  one  another.  For  some  events,  we  can  also  clearly 
identify  the  PcP  arrival. 

To  accurately  measure  the  differential  arrival  times  we  used  a  waveform  cross-correlation 
technique.  Cross-correlation  analysis  either  in  the  frequency  domain  (e.g.,  Poupinet  et  al.^ 
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Table  2:  Distances  and  Azimuths  Prom  JVE  Event  to  Stations 


Station 

Distance 

Azimuth  (deg. 

(deg) 

c.vj.  from  N) 

ANTO 

33.5 

-89.4 

BJI 

27.9 

96.2 

COL 

59.7 

20.9 

GBAR 

36.3 

-177.7 

GRFO 

42.1 

-63.1 

HIA 

26.2 

75.4 

KMI 

30.9 

134.4 

KONO 

39.2 

-48.4 

LZH 

22.7 

117.8 

YKAR 

67.1 

6.6 

1984,  Ito,  1985;  Fremont  and  Malone,  1987;  Moriya  et  al,  1994)  or  in  the  time  domain  (e.g., 
Prankel,  1982;  Pechmann  and  Kanamori,  1982;  Phillips  et  ai,  1992;  Deichmann  and  Garcia- 
Fernandez,  1992;  Rodi  et  ai,  1993;  Li  et  ai,  1995)  have  been  developed  to  quantitatively 
characterize  the  degree  of  similarity  of  seismic  waveforms  from  a  cluster  of  earthquakes  close  in 
space  and  to  measure  their  differential  arrival  times  in  an  accurate,  objective,  and  consistent 
manner.  The  conventional  time  domain  analysis  typically  enables  arrivals  times  to  be  read, 
at  best,  to  an  accuracy  of  one  sample  interval,  while  the  cross-spectral  method  (Poupinet  et 
ai,  1984,  Ito,  1985,  1990)  and  interpolation  techniques  in  the  time  domain  (Deichmann  and 
Garcia-Fernandez,  1992;  Li  et  ai,  1995)  can  improve  the  timing  precision  to  between  0,1  and 
0.5  sampling  intervals.  For  our  data  set  used  here,  the  accuracy  of  differential  arrival  time 
measurements  typically  ranges  from  0.01  to  0.02  seconds. 

Figure  6  is  an  example  of  the  waveform  cross-correlation  analysis  procedure.  The  top  frame 
of  Figure  6  shows  three  vertical  P  waveforms  of  three  Balapan  explosions  recorded  at  station 
BJI.  The  epicentral  distance  and  azimuth  are  28  and  96  degrees,  respectively.  The  waveforms 
are  aligned  by  the  origin  times  published  by  the  NEIC/USGS  (1987,  1988).  We  used  the  JVE 
explosion  (event  94)  as  a  master  event  and  explosions  85  and  86  (Table  1)  as  slave  events. 
The  cross-correlation  function  between  a  master  and  a  slave  event  was  calculated  for  a  P  wave 
window  to  measure  the  differential  arrival  time  between  events,  as  inferred  from  the  time  lag 
between  the  maximum  peaks  of  the  auto-  and  cross-correlation  functions.  The  P  wave  window 
length  used  here  is  6  s.  The  bottom  frame  of  Figure  6  shows  the  cross-correlation  functions  for 
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events  85  and  86,  as  well  as  the  auto-correlation  function  for  event  94.  The  peak  correlation 
values  are  0.95  and  0.9  for  events  85  and  86,  respectively,  suggesting  that  the  slave  events  (85, 
86)  are  indeed  located  very  close  to  the  master  event  94.  The  differential  arrival  times  relative 
.  to  event  94  are  0.475  and  -0.245  s  for  events  85  and  86,  respectively. 

Figure  7  shows  the  analogous  waveform  correlation  results  for  station  ANTO,  which  has  an 
epicentral  distance  of  34  degrees  and  an  azimuthal  of  271  degrees  from  the  test  site.  The  top 
frame  of  Figure  7  shows  the  vertical  component  P  waveforms  for  events  85,  86  and  94,  while 
the  cross-correlation  and  auto-correlation  functions  are  depicted  in  the  bottom  frame  of  Figure 
7.  The  P  wave  window  length  used  here  is  6  s.  The  cross-correlation  coefficients  are  again  very 
high,  0,92  to  0.94,  and  differential  P  wave  arrival  times  are  measured  to  be  about  -0.475  and 
0.36  s  for  events  85  and  86,  respectively,  relative  to  event  94. 

We  also  calculated  the  cross-correlation  functions  for  PcP  waves  at  some  stations.  Three 
vertical  component  seismograms  showing  both  P  and  PcP  waves  recorded  at  station  ANTO 
(distance  =  34  degrees,  azimuth  =  271  degrees)  are  shown  on  the  top  frame  of  Figure  8.  Using 
a  6  second  window  around  the  PcP  wave,  we  calculated  the  cross-correlation  functions  for  events 
85  and  86  and  the  auto-correlation  function  for  event  94  (bottom  frame  of  Figure  8).  The  cross¬ 
correlation  coefficients  are  0.8  and  0.72  for  events  85  and  86,  respectively.  Using  event  94  as  a 
reference,  the  PcP  differential  arrival  times  are  measured  to  be  -0.175  s  and  0.165  s  for  events 
85  and  86,  respectively. 

Table  3  lists  all  the  differential  times  that  were  determined  together  with  their  estimated 
accuracies  (a).  Figure  9  summarizes  the  results  of  the  waveform  correlation  analysis  (Figures 
6  to  8).  The  top  and  middle  frames  of  Figure  9  are  results  for  P  waves.  Two  stations  BJI  and 
ANTO  have  a  similar  epicentral  distance,  but  with  an  azimuth  difference  of  about  174  degrees 
they  are  in  opposite  directions  from  the  events.  For  event  85,  the  P  wave  differential  time  is 
0.475  s  at  station  BJI  and  -0.475  s  at  station  ANTO,  indicating  the  event  is  closer  than  the 
master  (94)  to  station  ANTO  and  farther  away  from  station  BJI.  This  observation  agrees  with 
the  determination  by  Thurber  et  a7(1993),  based  on  satellite  imaging,  that  event  85  is  located 
about  5  km  west  of  event  94  (Figure  2).  In  contrast,  for  event  86,  the  P  wave  differential  time 
is  -0.245  s  at  station  BJI  and  0.36  s  at  station  ANTO,  indicating  that  this  event  is  east  of  event 
94,  which  is  also  consistent  with  Thurber  et  al.  (1993). 
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Table  3;  Differential  Arrival  Times  Relative  to  JVE  Event 


Station 

Phase 

Slave  Event 

Time  (s) 

(7  (s) 

ANTO 

P 

85 

-0.475 

0.025 

ANTO 

PcP 

85 

-0.175 

0.03 

ANTO 

P 

86 

0.360 

0.01 

ANTO 

PcP 

86 

0.165 

0.015 

BJI 

P 

85 

0.475 

0.01 

BJI 

P 

86 

-0.245 

0.01 

BJI 

P 

88 

0.290 

0.01 

BJI 

P 

89 

0.765 

0.015 

COL 

P 

85 

-0.365 

0.015 

COL 

PcP 

85 

-0.360 

0.01 

COL 

P 

86 

-0.250 

0.01 

COL 

PcP 

86 

-0.215 

0.015 

COL 

P 

88 

-0.565 

0.015 

COL 

PcP 

88 

-0.520 

0.02 

GEAR 

p 

85 

0.660 

0.01 

GEAR 

PcP 

85 

0.340 

0.01 

GEAR 

p 

86 

0.315 

0.015 

GEAR 

PcP 

86 

0.315 

0.015 

GEAR 

p 

88 

1.000 

0.015 

GEAR 

PcP 

88 

0.480 

0.02 

GEAR 

p 

89 

0.120 

0.02 

GEAR 

PcP 

89 

0.060 

0.01 

GRFO 

p 

85 

-0.680 

0.02 

GRFO 

p 

86 

0.330 

0.02 

GRFO 

p 

88 

-0.410 

0.01 

GRFO 

p 

89 

-0.345 

0.01 

HIA 

p 

85 

0.262 

0.02 

HIA 

PcP 

85 

0.075 

0.01 

HIA 

p 

86 

-0.275 

0.015 

HIA 

PcP 

86 

-0.050 

0.01 

HIA 

p 

88 

0.063 

0.015 

HIA 

PcP 

88 

0.060 

0.01 

HIA 

p 

89 

0.775 

0.015 

HIA 

PcP 

89 

0.350 

0.01 

KMI 

p 

85 

0.600 

0.01 

KMI 

p 

86 

-0.112 

0.015 

KMI 

p 

88 

0.690 

0.01 

KMI 

p 

89 

0.425 

0.05 

KONO 

p 

85 

-0.770 

0.02 

KONO 

p 

86 

0.270 

0.02 

KONO 

p 

88 

-0.515 

0.015 

LZH 

p 

85 

0.685 

0.015 

LZH 

p 

86 

-0.260 

0.01 

LZH 

p 

88 

0.650 

0.01 

LZH 

p 

89 

0.715 

0.01 

YKAR 

p 

89 

0.360 

0.01 

YKAR 

PcP 

89 

0.285 

0.015 
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Table  4:  Locations  of  Five  Presumed  Underground  Nuclear  Explosions  Determined  from  Dif¬ 
ferential  Arrival  Times 


Event 

Latitude 

Longitude 

North  (km) 

East  (km) 

85 

49.9574 

78.7591 

8.8 

-4.6 

86 

49.9050 

78.8874 

3.0 

4.5 

88 

49.9786 

78.8052 

11.2 

-1.3 

89 

49.8601 

78.7117 

-2.0 

-8.0 

94 

49.8781 

78.8239 

0.0 

0.0 

Table  5:  90%  Confidence  Ellipses  on  Locations  Relative  to  Event  94 


Strike  (deg 

Semi-major 

Semi-minor 

Event 

c.w.  from  N) 

axis  (km) 

axis  (km) 

85 

42 

1.3 

1.1 

86 

18 

1.4 

0.9 

88 

31 

1.6 

0.9 

89 

30 

1.7 

1.0 

RELOCATION  RESULTS  AND  CONCLUSIONS 

We  applied  our  multiple  event  location  algorithm  (Rodi  et  a/.,  1994)  to  the  differential  arrival 
time  data  listed  in  Table  3.  The  epicenters  and  origin  times  of  the  four  slave  events  (85,  86, 
88,  89)  were  fit  to  the  data  in  a  weighted  least  squares  sense.  The  epicenter  and  origin  time  of 
event  94  were  fixed  to  the  values  given  in  Table  1  (Thurber  et  aZ.,  1993;  NEIC/USGS,  19987, 
1988).  The  depths  of  all  events  were  constrained  to  the  earth’s  surface.  The  forward  model  for 
traveltimes  was  computed  from  the  IASP91  tables,  retrieved  from  the  Center  for  Monitoring 
Research  in  Arlington,  Virginia.  Our  relocated  epicenters  are  listed  in  Table  4  and  plotted  in 
Figure  10. 

The  final  residuals  of  fit  to  the  differential  arrival  time  data  averaged  5.4  standard  deviations. 
Allowing  for  degrees  of  freedom,  this  implies  posterior  estimates  of  the  data  standard  deviations 
equal  to  6.2  times  the  assumed  (prior)  standard  deviations  shown  in  Table  3.  We  attribute  the 
large  residuals  to  two  major  causes.  First,  the  events  were  fixed  to  a  common  depth  while 
the  true  depths  of  the  events  vary.  Second,  the  IASP91  traveltime  tables  do  not  reflect  the 
relatively  lower  velocity  layers  at  shallower  depths  at  the  KTS  site  (Priestley  et  a/.,  1988;  Li 


33 


! 

i 

and  Thurber,  1991;  Quin  and  Tliurber,  1992);  nor  do  they  reflect  lateral  velocity  variations  at 
the  site. 

Figure  10  compares  our  locations  (stars)  with  the  NEIC/USGS  locations  (octagons)  and 
satellite  locations  of  Thurber  et  al  (1993)  (squares).  The  USGS  locations  have  been  shifted 
so  as  to  align  event  94  with  its  satellite  location.  Our  location  for  event  94  coincides  with  the 
satellite  location  as  a  constraint.  We  see  that  our  locations  for  the  remaining  four  events  are 
within  2  to  3  km  of  the  satellite  locations  of  Thurber  et  al,  (1993).  Our  location  for  event  85  is 
comparable  to  the  NEIC  location  for  this  event,  but  for  the  other  three  events  (86,  88  and  89) 
our  locations  are  better  than  the  NEIC  locations,  with  mislocations  reduced  by  factors  of  about 
2  to  5.  We  also  note  that  if  we  had  used  more  accurate  traveltime  tables,  allowing  for  slower 
shallow  velocities,  our  pattern  of  locations  would  shrink  and  its  agreement  with  the  satellite 
image  locations  would  improve  slightly.  In  Table  5  we  list  the  formal  confidence  ellipses  for 
the  location  of  each  slave  event  relative  to  the  location  of  the  JVE  event  (94).  (These  error 
estimates  include  the  posterior  variance  factor.)  We  see  that  the  90%  confidence  ellipses  have 
semi-axes  of  about  1  and  1.5  km,  which  is  about  half  of  the  actual  location  errors. 

We  conclude  from  these  results  that  a  relative  event  location  method  applied  to  differential 
P  and  PcP  arrival  times  at  a  small  number  of  far-regional  and  teleseismic  stations  yields  relative 
location  accuracies  for  explosions  as  good  as,  or  better,  than  conventional  locations  obtained 
with  hundrededs  of  arrival  time  picks  from  a  large  network  of  globally  distributed  stations. 
The  reason  for  this  is  that  the  waveforms  of  explosions  differing  in  location  by  only  several 
kilometers  are  very  similar,  allowing  the  measurement  of  highly  accurate  arrival  time  differences 
between  events  using  cross-correlation  techniques.  Therefore,  waveform  analysis  and  relative 
event  location  methods  can  be  valuable  assets  for  nuclear  monitoring  and  discrimination. 
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Figure  1:  World  map  showing  the  locations  of  five  presumed  nuclear  explosions  (octagons)  in 
Balapan,  Kazakhstan  test  site  (KTS)  and  ten  seismic  stations  (triangles)  with  epicentral 
distances  ranging  from  23  to  67  degrees. 
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Figure  2:  Map  of  Kazakhstan  test  site  (KTS)  showing  relative  locations  of  five  presumed  nuclear 
explosions  determined  by  Thurber,  Quin  and  Richards  (1993)  using  a  satellite  imaging 
technique  (squares),  and  by  the  USGS  with  hundreds  of  arrival  time  picks  (octagons). 
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Bolopon  Explosion  94  (880914-03:59:57.4)  Recorded  at  10  Stations 


Time  (s) 


Figure  3:  Short-period,  vertical  component  seismograms  of  the  14  September  1988  Joint  Veri¬ 
fication  Experiment  (JVE)  explosion  (event  94)  at  KTS  recorded  at  ten  seismic  stations  at 
far-regional  and  teleseismic  distances. 
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Balapan  Explosions  Recorded  at  Station  LZH  (D=23,  AZ=118) 
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Figure  4;  Vertical  component  seismograms  of  five  presumed  nuclear  explosions  recorded  by 
station  LZH  at  a  far-regional  distance.  Distance  (D)  and  azimuth  (AZ)  are  given  in  the 
title  in  degrees. 
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Balapan  Explosions  Recorded  at  Station  COL  (D=60,  AZ=21) 
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Figure  5:  Vertical  component  seismograms  of  five  presumed  nuclear  explosions  recorded  by 
station  COL  at  a  teleseismic  distance. 


Figure  6:  Top  frame:  P  waves  on  vertical  component  seismograms  of  three  explosions  at  KTS 
recorded  at  station  B  JL  Bottom  frame:  P  wave  cross-correlation  functions  for  events  85  and 
86  and  auto-correlation  function  for  event  94  (master  event). 
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Bolopan  Explosions  Recorded  ot  Station  BJI  (D=28,  AZ=96) 


Figure  7:  Top  frame:  P  waves  on  vertical  component  seismograms  of  three  explosions  at  KTS 
recorded  at  station  ANTO.  Bottom  frame:  P  wave  cross-correlation  functions  for  events  85 
and  86  and  auto- correlation  function  for  event  94  (master  event). 


Balopan  Explosions  Recorded  at  Station  ANTO  (D=34,  AZ=270) 


PcP  Wove  Cross-Correlation  Analysis:  ANTO  (D=34,  AZ=270) 


Figure  8:  Top  frame:  P  and  PcP  waves  on  vertical  component  seismograms  of  three  explosions 
at  KTS  recorded  at  station  ANTO.  Bottom  frame:  PcP  wave  cross-correlation  functions  for 
events  85  and  86  and  the  auto-correlation  function  for  event  94  (master  event). 
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Figure  9:  Comparison  of  cross-correlation  functions  at  different  stations  and  calculated  with 
different  phases.  Top  frame:  P  wave  cross-correlation  functions  calculated  for  events  85  and 
86  at  station  BJI.  Middle  frame:  P  wave  cross-correlation  functions  calculated  for  events 
85  and  86  at  station  ANTO.  Note  the  diflFerences  between  the  P  wave  differential  arrival 
times  for  the  same  event  at  different  stations.  Bottom  frame:  PcP  wave  cross-correlation 
functions  calculated  for  events  85  and  86  at  station  ANTO. 
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Figure  10:  Comparison  of  our  relocation  results  (stars)  with  the  locations  determined  with 
satellite  image  analysis  (squares)  by  Thurber,  Quin  and  Richards  (1993),  and  with  hundreds 
of  absolute  arrival  time  picks  by  USGS  (octagons). 
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ABSTRACT 


new  method,  which  we  name  Spherical  Mapping  Approximation  (SMA),  is  de- 
ve  oped  for  the  evaluation  of  displacement  fields  of  body  waves  and  surface  waves 
rom  explosions  in  non-spherical  cavities  embedded  in  elastic  media. 


nder  SMA,  the  explosion-generated  stress  distribution  on  the  surface  of  an  arbi¬ 
trary  cavity  is  mapped  onto  the  surface  of  an  equivalent  virtual  spherical  cavity  having 
t  e  volume  of  the  true  cavity.  The  analytical  results  express  the  displacement  field  in 
terms  of  a  multipole  double-sum  expansion  of  spherical  eigenvectors  with  coefficients 
in  the  form  of  a  finite  Legendre  transform  of  the  components  of  the  normal  vector  of 
the  cavity  boundary.  These  ’cavity  integrals’  can  be  evaluated  exactly  for  spheroidal 
and  cyllindrical  inclusions. 

In  the  long- wave  far-field  approximation,  symmetric  finite  cavities  are  shown  to 
e  equivalent  to  a  linear  combination  of  point  dipoles  directed  along  the  principal 
cavity-axes.  The  ensuing  radiation  patterns  yield,  in  general,  4-lobe  patterns  for  S- 
waves,  two-lobe  patterns  for  P-waves  and  single  to  two-lobe  Rayleigh-wave  patterns, 
independent  of  the  details  of  the  cavities’  shape.  However,  all  radiation  patterns 
are  modulated  by  a  frequency-dependent  ’cavity-factor’  that  embodies  the  boundary 
conditions  on  the  cavity  surface. 

Moreover,  it  is  shown  that  the  radiation  pattern  for  P-waves  from  a  non-spherical 
symmetrical  cavity  in  the  long-wave  far-field  approximation  is  always  dipolar.  Since 
the  radiation  pattern  of  radiated  P-waves  from  a  standard  earthquake  is  always 
quadrupolar,  the  cavity  explosion  behaves  like  a  non  double-couple  earthquake.  Thus, 
the  examination  of  the  deviatoric  moment  tensor  of  a  given  seismic  event  enables  one, 
in  principle,  to  state  whether  it  is  a  standard  earthquake  or  perhaps  (if  the  S-wave 
pattern  is  quadrupolar)  an  evasion  of  the  test-ban  treaty. 

Displacement  patterns  for  body  and  surface  waves  are  calculated  for  spheroidal 
an  cylindrical  cavities  for  a  wide  range  of  aspect  ratios  and  corresponding  aperture 
angles,  exhibiting  the  whole  range  of  cavity  shapes  from  a  line  source  to  a  disc.  The 
moments  of  the  equivalent  dipoles  are  shown  to  depend  on  the  corresponding  cavity- 
integrals,  the  elastic  constants  of  the  medium  in  the  neighborhood  of  the  source,  and 
the  initial  energy  injection. 

All  non-spherical  cavities  generate  strong  shear-waves,  except  for  special  aperture 
angles  at  which  a  spherical  P-wave  is  generated,  unaccompanied  by  S-waves. 

The  wave-spectra  of  body  waves  (surface  waves)  exhibit  a  corner-frequency  (peak 
frequency)  at  a  wavelength  equal  to  the  radius  of  the  equivalent  sphere.  This  enables 
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one  to  deduce  the  size  of  the  cavity  from  the  spectrum  of  its  far-field  displacement 
signals,  provided  that  the  explosion  is  fully  decoupled  and  that  the  interaction  of  the 
shock  wave  with  the  medium  occurred  in  the  elastic  regime. 

The  results  of  the  present  research  are  applicable  to  the  detection  and  identifica¬ 
tion  of  seismic  signals  from  clandestine  underground  nuclear  explosions. 


1.  Introduction 


The  advent  of  testing  underground  nuclear  explosions  during  1957-1966  gave 
rise  to  intensive  research  efforts  on  the  subject  of  seismic  wave  radiation  from  ex¬ 
plosions  in  spherical  cavities.  The  concept  that  it  might  be  possible  to  significantly 
reduce  the  radiated  seismic  signal  by  detonating  the  device  in  a  cavity  (known  as 
decoupling),  was  first  proposed  publicly  by  A.L.  Latter  at  the  1959  Nuclear  Test 
Ban  Conference  in  Geneva  {Latter  et  aL,  1961).  The  possibility  of  concealing  nu- 
c  ear  explosions  then  became  a  vital  issue  in  connection  with  the  test-ban  treaty, 
and  its  importance  has  presently  increased  due  to  the  new  treaty  to  be  signed  and 
ratified  internationally  in  1996. 

In  spite  of  the  existence  of  a  vast  literature  of  the  subject  of  decoupling  [Lat¬ 
ter  ei  al.,  1959;  Latter  et  ai.,  1961a, b;  Haskell,  1961;  Herbst  et  ai,  1961;  Patter- 
Springer,  1966;  Werth  and  Randolph,  1966;  Rogers,  1966;  Lewin  and 
Trehnan,  1966;  Rawson  et  d.,  1966;  Haskell,  1967;  Springer  et  a/.,  1968;  King  et 
d.,  1989;  Stevens  et  d.,  1991;  Adushkin  et  d.,  1993;  Florence  and  Miller,  1993; 

Glenn,  1993),  a  number  of  major  issues  of  importance  w.r.t.  sei.smic  monitoring 
still  remain  unsolved. 

Elastodynamical  problems  associated  with  spheres  and  spherical  cavities  were 
^Ived  throughout  the  Ip''*  century  by  many  mathematical  physicists  (Love,  1927). 

owever,  the  first  application  of  a  problem  of  this  type  to  seismology  was  made 
o^y  m  1942  by  Sharpe;  he  gave  explicit  expressions  for  the  time-domain  ground 
splacem^ents  produced  by  chemical  explosion  pressure  in  a  finite  cavity  embedded 
in  an  unbounded  elastic  solid. 

Sharpe  s  solution  assumes  total  spherical  symmetry  of  the  source  and  the 

medium.  It  therefore  predicts  a  radial  displacement  everywhere  and  the  absence 
of  shear  waves. 

However,  seismograms  produced  by  underground  nuclear  explosions  are  con¬ 
siderably  more  complicated  than  what  one  would  expect  from  Sharpe’s  simple 
model  (e.g.  Toksoz  et  ai,  1964;  Johnson,  1988;  Taylor  et  ai,  1991). 

Numerous  attempts  were  made  during  the  past  four  decades  to  account  for 
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these  observations  in  terms  of  media  complexities  (iahomogeneity,  anisotropy,  pre¬ 
stress.  nonlinearity)  or  source  complexities  (asphericity,  spall).  But  a  serious  at¬ 
tempts  to  extend  Sharpe’s  analytical  solution  beyond  that  of  a  spherical  cavity  was 
lacking.  For  this  rerison,  all  scaling  laws  used  to  date  (King  et  a/.,  1989;  Stevens  et 
a/.,  1991;  Adushkin  et  a/.,  1993;  Florence  and  Miller,  1993)  to  estimate  the  yield 
of  nuclear  explosions  from  the  seismic  data  are  still  based  on  a  modification  of 
Sharpe’s  solution  given  by  Haskell  (1961,  1967). 

In  surveying  the  literature  on  the  subject  since  1942,  we  could  discern  two 
main  avenues  of  approach: 

I.  Helmholtz-type  integral  formulas:  These  illustrate  Huygen’s  principle  for 
the  two  wave  fronts  of  the  elastic  wave  field  via  a  vector  Green’s  theorem.  It 
is  known  also  as  the  ‘Representation  Theorem’  (Ben-Menahem  and  Singh,  1981, 
p.  174)  and  states  that  the  displacement  field  u  outside  the  source  region  is  given 
by  two  surface  integrals  over  the  confining  cavity  5 

a(r)  =  J  F(0  ■  G(r  I  ()dS(0  -  J  u(0n(0:  T^{G(r  \  0]dSi0.  (1) 

s  s 

Here  ^  are  source- coordinates,  n(^)  is  the  outward  unit  noi'inal  vector  to  S,  F{^)  = 

n(()  •  T(u)  is  the  normal  traction  vector,  dS  is  the  surface  area  element,  G'(r  |  ^) 

IS  the  Green’s  tensor  for  the  elastic  medium  (finite  or  infinite),  T (u)  =  A  /  div't7  + 

3 

-f  uV)  is  the  stress  dyadic,  while  T  is  a  stress  tryadic  (third  rank  Green’s 
stress  tensor).  It  is  given  explicitly  in  either  index-free  or  indicia!  form 

T^(r  I  0  =  A?  div^  G(r  |  f)  +  m(  G  -f  G^  V), 

3 

mri  ■)"  dm  Gtn). 

where  dm  = 

Since  F(^)  is  usually  given  and  G  is  known,  tlie  first  integral  of  (1)  can  be 
evaluated,  analytically  or  numerically,  on  S.  The  second  integral,  requires  the 
preknowledge  of  the  displacements  on  S.  Since  in  most  cases,  this  function  is 
apriori  unknown,  it  must  be  first  calculated  from  (1)  itself,  acting  to  this  end 
as  an  integral  equation  for  u  on  S.  This  drawback,  apart  from  the  nontrivial 
integration  over  S,  makes  this  approach  quite  tedious. 

Varadan  (1976),  mitigated  the  hardships  of  the  integral-formula  approach  by 
applying  to  it  the  method  of  eigenfunction-expansion,  well-known  from  the  theory 

of  integral  equations  (e.g.  Kondon,  1991).  He  expands  G,  ^xnd  T{G)  in 
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terms  of  a  spherical  eigenvector  base,  and  solves  for  the  unknown  coefficients  of 
the  expansion  of  the  displacements  u(f).  An  infinite  matrix  is  involved,  but  in 
many  cases  good  accuracy  can  be  obtained  with  a  truncated  matrix  of  60  x  60. 

Glenn  et  al.  (1985)  and  Rial  and  Moran  (1986)  solved  (1)  numerically  for  the 
surface  displacement  over  a  spheroid.  A  finite  element  code  was  used  to  compute 
the  displacement  on  the  cavity  surface  and  then  (1)  was  used  again  to  find  the 
displacements  in  the  far-field.  No  analytic  solution  was  given.  But  the  singular 
nature  of  the  integral,  when  the  observer  point  is  on  the  surface,  makes  it  difficult  to 
apply  this  scheme  to  arbitrary  cavity  geometries.  Also,  the  finite-element  solution 
is  computationally  inefficient  for  direct  calculation  of  the  far-field. 

Stevens  et  al.  (1991)  reported  non-analytic  finite  difference  simulations  of 
partially-coupled  explosions  in  an  ellipsoidal  cavity  having  an  aspect  ratio  of  4  to 

I.  Since  their  report  includes  no  visible  mathematics,  except  for  a  fiat  statement 
that  they  have  integrated  Eq.  (1)  numerically  in  the  time  domain  for  the  far-field, 
a  straightforward  evaluation  of  their  results  is  not  fcasil:>le. 

II.  Boundary-value  solutions  for  non-spherical  configuration:  ffeeian 
(1953)  has  tried  to  solve  the  problem  of  radiation  from  a  cylindrical  cavity  in  an 
isotropic  infinite  elastic  space.  He  assumed  a  cylindrical  hole  of  infinite  length  with 
a  prescribed  stress  applied  to  a  finite  length  of  its  wall,  which  is  not  a  cavity.  But 
even  so,  his  solution  is  plagued  with  serious  analytical  inaccuracies  (Aljo-Zena, 
1977;  Usami  and  Hirono,  1956).  Hazebroek  (1966)  considered  a  special  case  of 
a  finite  line  source  given  as  a  limit  of  a  narrow  elongated  ellipsoid  of  revolution 
whose  minor  axis  tends  to  zero.  The  interior  surface  of  this  ellipsoid  is  subjected 
to  a  pressure  which  is  a  given  function  of  the  time  only.  It  was  concluded  that  the 
shape  and  magnitude  of  the  compressional  waves  were  independent  of  direction 
and  that  the  magnitude  of  the  shear  wave  is  maximum  at  an  angle  of  45°  to  the 
direction  of  the  line  source. 

Usami  and  Hirono  (1956)  and  Usami  (1958)  considered  the  elastic  wave  gen¬ 
erated  from  prolate  and  oblate  spheroidal  cavities  whose  walls  were  subjected  to 
normal  harmonic  stress.  The  above  authors  used  si)heroidaI  coordinates  from  the 
start,  and  their  solution  is  therefore  expressed  in  terms  of  the  sphei'oidal  Baer 
eigenfunctions  (Moon  and  Spencer,  1971)  which  are  eventually  expanded  in  terms 
of  spherical  Bessel  functions.  Their  numerical  calculations  revealed  for  the  ratio 
d  —  {  sou^c^rre^ngth^^^^  }  ~  3,  a  4-lobe  radiation  pattern  for  compressional  waves 
in  the  far-field. 

Zhao  and  Harkrider  (1992)  gave  an  analytic  formulation  for  the  wave  fields 
from  an  off-center  explosion  in  an  embedded  solid  sphere  in  an  elastic  whole-space. 
Their  calculations  show  that  the  degree  of  shear-wave  generation  is  determined  by 
the  asymmetry  of  the  source  region.  The  radiation  patterns  at  different  periods 
for  different  parameters  of  the  media  suggest  that  the  asymmetry  of  the  source 
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region  has  significant  effects  on  spectral  components  with  the  ratio  d  <1.  Their 
model  was  intended  to  simulate  a  tamped  nuclear  explosion. 

A  more  realistic  model  of  an  underground  nuclear  explosion  in  a  preexisting 
cavity  requires  the  presence  of  a  stress-free  planar  boundary  to  accommodate  for 
the  earth’s  surface.  Analytically,  this  problem  is  of  a  higher  degree  of  complexity 
since  it  involves  the  simultaneous  use  of  two  different  coordinate  system  —  the 
spherical  and  the  cylindrical. 

The  first  attempt  to  solve  this  difficult  problem  is  due  to  Dcn-Menahem  and 
Cisternas  (1963),  who  used  the  ‘Erdelyi  Integral  Formula’  to  transform  the  elastic 
fields  between  the  two  said  coordinates.  The  convergence  of  their  solution  wtis  later 
established  by  Thiruvenkatachar  and  Viswanathan  (19G5)  and  Gregory  (1967). 

The  present  paper  offers  a  new  simplified  approach  which  is  especially  suit¬ 
able  to  the  physical  conditions  prevailing  in  underground  nuclear  explosions.  The 
basic  assumption  (supported  by  observations  over  the  past  forty  years)  states  that 
irrespective  of  the  shape  of  the  cavity  or  the  point  at  which  explosive  device  is 
placed  in  it  —  the  direction  of  pressure  on  its  walls  is  always  normal  to  the  walls 
at  every  point,  such  that  there  is  no  initial  tangential  force.  The  pressure  on  the 
walls  is  activated  practically  simultaneously  at  all  points  of  the  boundary. 

We  analyze  comparatively  three  typical  situations,  listed  in  increasing  order 
of  complexity  (Fig.  1): 

(a)  Source  at  the  center  of  a  spherical  cavity.  The  initial  traction  vector  on  the 
wall  is  radial  and  uniform  at  all  points.  No  shear  waves  arc  produced  in  the 
outside  homogeneous  elastic  medium  at  any  time  (Fig.  la). 

(b)  Source  is  at  the  center  of  a  spheroidal  cavity  (Fig.  lb).  Since  the  normal 
is  not  in  the  radial  direction  from  the  center,  shear  -waves  are  produced  ab 
initio  in  the  elastic  medium.  One  can  then  simulate  the  problem  in  terms  of 
a  mapping  of  the  stress  distribution  on  the  surface  of  an  cquiualcnt  virtual 
spherical  cavity.  Details  of  this  process  are  given  in  Section  2. 

(c)  A  stress-free  planar  boundary  is  introduced  to  model  an  underground  explo¬ 
sion  in  a  tunnel  whose  major  axis  is  parallel  to  the  free  surface  (Fig.  Ic). 

With  this  step  accomplished,  the  source-field  of  an  oblate  spheroid  in  an  in¬ 
finite  medium  is  obtained^.  The  second  step  included  the  integral  transformation 
of  the  spherical  eigenvectors  into  cylindrical  eigen^x^ctors,  appropriate  for  the  new 
boundary  conditions  needed  over  the  planar  boundal•J^  Once  the  l)oundary  condi¬ 
tions  are  stated,  the  residues  at  the  poles  of  the  integral  expressions  are  evaluated 
in  order  to  obtain  a  closed-form  expressions  for  the  E.ayleigh  wave  field. 

The  results  of  the  above  cases  are  calculated  and  exhibited  in  graphical  form 

^  By  ‘oblate^  we  mean  here  that  the  cavity  has  a  horizont  al  axis  of  syinnietry. 
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which  emphcLsize  the  dependence  of  the  field  on  the  explosion  geometry  and  source 
frequency.  The  shape  of  the  cavity  strongly  affect  the  ensuing  radiation  patterns. 

The  new  idea  of  the  equivalent  spherical  ca^'ity  is  very  helpful  in  bringing 
about  the  essential  features  of  the  resulting  elastodynamic  field  that  is  transmitted 
outwards  from  the  cavity.  Clearly,  we  are  not  able  to  quantify  the  entire  complex 
physical  process  of  an  underground  nuclear  explosion,  due  to  lack  of  measurements 
in  the  cavity  and  its  adjacent  neighborhood  during  the  first  few  milliseconds  or 
so  after  the  explosion.  But  practically  we  are  interested  onl}"  in  that  tiny  fraction 
of  energy  that  is  converted  into  seismic  waves.  Since  the  boundary  conditions 
associated  with  this  conversion  are  only  vaguely  known,  there  is  no  sense  in  being 
too  exact  in  the  mathematical  elastic  model.  Thus,  instead  of  giving  an  exact 
mathematical  solution  to  an  approximate  physical  situation,  we  prefer  to  give  an 
approximate  mathematical  solution  to  an  equivalent  physical  model,  which  we 
have  precisely  formulated. 


2.  The  Spherical  Mapping  Approximation 

Consider  first  a  spherical  cavity  of  radius  r  =  vq  that  is  subjected  to  a  pre¬ 
scribed  distribution  of  stresses  on  its  inner  wall.  To  estal:)lish  the  displacement 
field  u  at  an  observation  point  r  outside  the  cavity  we  assume  an  expansion  of  the 
field  in  terms  of  the  Hansen  spherical  eigenvectors  (Appendix  A), 

,  (3) 

where  kp  =  ks  ^  and  u  is  the  angular  frequency.  Let  the  boundary 
conditions  at  the  surface  of  the  cavity  be 

er  •  T{u))  —  F{6,(p)  at  r  =  tq,  (4) 

where  is  a  unit  vector  in  the  direction  of  increasing  r,  and  T (it)  =  AJ divu  + 
/i(Vu  +  uV)  is  the  stress  dyadic  in  a  medium  ha\’iiig  the  Lame  coefficients  A,  fi. 
Since  .F  is  known,  it  can  be  expanded  into  a  series  of  vector  spherical  harmonics 
{rme^Bm£,Cme}  with  determinable  coefficients  )  ( Ben-Menahem 

and  Singh,  1981,  pp.  221-222) 

(5) 

m,£ 

The  explicit  expressions  for  the  vectors  ^rnC^^mc}  is 

given  in  Appendix  A. 
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Inserting  u  from  (3)  into  (4),  straightforward  calculations  lead  to  the  explicit 
expressions  for  the  unknown  displacement  coefficients  in  terms  of  the  known  stress 
distribution 


Pml  = 


fiks 


1  -^L-^f,2(x)  +  2^(^+l)7LFf,](x) 

2nkp  Ae 


(6) 

(7) 


where 


X  =  k^ro,  C  = 

=  2£(<!+  l)F^,i(x)F^,i(C)  - 


f<,2(x)  = 


Fi,3{^)  = 


:{i^  -  1)  -  1 


2i 


h^f  H-i  )  H — k\’.{x). 


2  ,.(2) 


f+1' 


=  spherical  Hankel  function  of  the  second  kind.  J 


(S) 


The  above  analysis  can  be  applied  to  model  an  explosion  in  a  preexisting  cavity; 
whatever  the  geometrical  shape  of  the  cavity,  the  explosion  gas  is  assumed  to 
exert  on  its  walls  a  uniform  pressure  equal  to  pQ{uj)n,  where  n  is  a  unit  normal 
to  any  given  point  and  po(^)  is  the  Fourier  transform  of  the  time- dependence  of 
the  initial  stress-pulse.  In  the  case  of  a  sphere  il  =  e,.  is  the  unit  radial  vector 
in  spherical  coordinates.  In  the  general  case  n  =  V5/1VS|,  where  S{x,y,z)  =  0 
is  the  equation  of  the  cavities’  surface.  Because  the  Ijoundary  conditions  require 
that  the  normal  stress  Fn  at  the  surface  of  the  cavity  be  etpial  to  the  negative  of 
the  applied  pressure,  we  write 


Fn  [r(0,  if)-,  6,  ip]  =  -po{uj)n 

=  -po(a')  [/(^,  (p)er  +  (j{O,ip)c0  +  h{6,<p)e^], 

where  {r,9,ip]  are  the  spherical  coordinates  of  a  point  on  S  (Fig.  1),  and  {/,</,  h} 
depend  on  the  geometry  of  the  surface  S.  Since  the;  surface  of  the  cavity  after  the 
explosion  does  not  remain  intact  due  to  non-elastic  deformations,  the  physics  of 
the  problem  can  accommodate  a  mathematical  simplification:  instead  of  applying 
the  boundary  conditions  on  the  aspherical  cavity,  we  map  tlu'iii  on  the  walls  of 
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a  virtual  spherical  cavity  which  has  the  same  voluTnc  as  the  real  cavity.  The 
geometry  of  the  original  cavity  still  enters  via  the  angular  dependence  of  n. 

With  this  provision,  the  mathematical  apparatus  expounded  in  (l)-(7)  can 
be  applied  directly.  Since  in  nuclear  explosions  the  distances  involved  in  changing 
from  S  to  the  sphere  r  =  tq  are  of  the  order  of  a  few  meters  only,  the  travel-time 
error  incurred  will  be  of  the  order  of  a  few  milliseconds,  which  is  totally  negligible. 
Hence,  we  are  able  to  avoid  the  use  of  horrendous  integrations  over  aspherical 
surfaces  as  required  by  other  methods  (e.g.  Varatharajulu  and  Pao,  1976). 

Let  5(r,  0,  ip-,e)=0  represent  a  smooth  surface  in  three  diincnsional  Euclidean 
space,  where  (r,  9,  ip)  are  the  spherical  coordinates  of  a  general  point  of  S  relative  to 
some  origin  (usually  chosen  inside  5)  and  0  <  e  <  1  is  a  dimensionless  geometrical 
pcirameter.  We  shall  assume  that  the  components  of  V5,  together  with  their  first 
and  second  derivatives  are  continuous  functions  of  9  and  cp  in  the  range  0  <  9  <  tt, 
0  <  (p  <  2x. 

The  explicit  expression  for  the  normal  to  5  is 


where 


V5 

|V5| 


=  +  g{9,p;e}c()  +  h{9,p: 


_J_^  1  dS 

A  dr '  ^  Ar  89 '  *  ArsinSd-p' 


\  /  r^  \d9  )  r^  sin*^  9  \  dp> 


2-,  1/2 


(10) 


(11) 


For  a  given  surface  5,  such  that  {fiQih)  together  with  their  first  and  second 
derivatives  are  continuous  functions  of  9  and  ^p  over  the  entire  range  of  (^,  </’),  we 
may  seek  an  expansion  of  the  vector  normal  n  in  terms  of  tin;  vector  spherical 
harmonics  (Appendix  A) 


oo  I 


"  =  E  E  [fmi  Pmt  +  9mi  \/ +1)  PmC  +  ^hncV +  1 )  C\nc]  • 


(12) 


Using  the  orthogonality  relations  of  the  scalar  spherical  hariuonics,  the  coefficients 
hmi}  are  explicitly  recoverable  in  terms  of  the  hiouni  functions  {/,</,  h} 


fmi  — 


2e  +  i{e-m)\ 


27r 


9mi  — 


47r  (C  +  m)! 

0  0 

2tt  tz 

2^+\  (£-m)! 


j  dtp  j  f(9,  tp-  e){  P/'‘(cos  9)t-^"^'^  }  sin  9d9.  (13) 


47r^(^  -f  1)  (^  +  m) 


dp 


0  0 
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(14) 


X  {PP {cos sin 0d9, 


hmi  = 


2n  TV 

2e  +  l  {i~m)\  /■  ,  //  1  ,0  d 

0  0 

X  {Pp{cose)e-^'^'^}smdde. 


47r£(£  +  l)(^  +  m)!7  “^7  '^^-'^"'09 


(15) 


In  the  special  case  of  azimuthally-symmetric  surfaces  =  0,  h  =  O] ,  the  above 
results  simplify  to 


/£  = 


2i+l 


TV 

j  £)Fg(cos  0}smB(W, 


(16) 


9l  — 


2i+l 

2i{i+l) 


j  g{9,if,e)-^Pi{cos9)sin9d9,  ^’ >  1,  (17) 


90  =  0. 


Applying  the  stress  distribution  given  in  (9)-(10)  to  (C),  using  (C)  (S),  we  find  by 
comparison,  for  edl  values  of  i  and  m 


/Sm<  =  tI;  =  -M^)Sme,  «“.(  =  -P(l(‘^)*-..f  ■ 


(18) 


The  final  expressions  for  the  respective  compressiunal  (P)  and  shear  (5)  displace¬ 
ment  fields  at  locations  r  >  tq,  are  obtained  via  (3)-(8) 

-  ^  po(t^)  /m£-F£,2(^s’'o)  -  2^(^  +  l)</;»if-ff,l(^'tsn)) 


2  pc /up 


m.i 


u 


,{r-,uj)  =  - 


P0(w) 


y 

,[{^sro)Ffi{ksro) 

m,t 

,  fmiFt^likp^o)  -  9mlF'e,3{kprQ)  ^ 


(20) 


For  azimuthally-symmetric  surfaces,  (19)-(20)  degenerate  into  the  simpler  form 

Po(‘^) 


Up{f,u)  = 


2fj,k 


Ug{f,u)  =  - 


E 

P  C=0,2A,... 

X  LQ^{kprQ), 

Po(^) 

1^’=’  (=(1,2,4,... 
X  N~{ksr), 


Fioiksf)  Fc  {{^'st) 

fd^)  ’  .  ■■-2C{C  +  l)gc{e)- 


E 


A( 


j.  f  ^F'i  i{kpr)  ^F^  -j{kpr) 

fd^) — X - 9d^) 


Af 


(21) 


A£ 


Af 
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(22) 


where  JVq  Q(^'s^)  =  0. 

The  above  equations  are  recastable  in  the  explicit  coniponcnt-fonn: 


P-waves: 


Ur  -  ^^^y]Ai{ksrQ-,e)Pi{cosB)Ii^'^\kpi-) 


UQ  = 


Po(^) 

fikp 


Y^At{ksro]e) 


dPlicosO)  li^c^\kpr) 


de 


kpV 


(23) 


S- waves: 


Xlr  — 


ug  = 


poM 

jikp 

Po(^) 

fikp 


'^Bi{kprQ;€)C{e  +  l)P£(cost^) 


hf\k,r) 

ksv 


'Y^Bi{kpro;e) 


dPc{cos  9) 


k 


■ar 


(24) 


where 


An  ^  ^Fe,2{^sro)  ,  F£.i(A-.s7-o) 


Bi{kprQ  \  e)  — 


2A£ 

r  ,  s.F't^iikpTQ)  Pc:dh''()) 

fii^) - T - gd^) 


A/ 


A, 


^'p 

k. 


(25) 

(26) 


For  a  spherical  cavity  gi  =  0,  fi  =  S^q  and  (21)-(22)  reduce  to  the  well  known 
result  (Ben-Menahem  and  Singh,  1981,  p.  222): 


_  P0{oj)  1  r-  n  .  ^ 

2Mip 

This  expression  can  also  be  recast  in  the  more  explicit  form 


(27) 


Up(f,u))  =  —  grad 


po(u7)a3  a..2c'-'/-o 


4P  + 


(28) 


=  —  grad 


where  ^'(t  — »•  oo)  =  ^(cj  — »■  0)  =  ,  ujq  =  ^  and  v{t  — +  oo)  =  (static 

deformation  field). 
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3.  Application  to  Realistic  Cavity  Shapes 


We  next  apply  our  ideas  to  the  simplest  departure  from  a  sphere  that  is 
adequate  to  model  underground  tunnels.  A  spheroidal  cavity  with  semi-major 
axis  c  (directed  along  the  z  axis)  and  semi-minor  axes  equal  to  a  is  surrounded 
by  an  infinite  homogeneous  and  isotropic  elastic  solid  with  Lame  coefficients  A, 
density  p  and  a  Poisson  ratio  a.  The  origin  of  the  coordinates  coincides  with  the 
spheroid’s  center  (Fig.  2). 

Let  the  spheroidal  surface  be  given  by  the  equation 

5’(i,y,2)  = -4(a:^ -f  y^)  + -4-“  -  1  =  0,  (29) 


where  c  >  a  is  the  major  semi-axis,  and  (x,y,z)  are  the  cartc'sian  coordinates  of 
any  point  P  on  the  surface.  The  normal  at  the  surface'  is  giveni  l)y  the  unit  vector 
n  =  V5/|VS'l.  Expressing  this  vector  in  spherical  coordinate's  at  P  (r,  we 

obtain,  after  simple  analysis 


n  =  erf{d;  e)  -h  e^y(^;  e);  e  =  -  <  1,  f  '^  +  (/^  =  1, 

c 

where  {er,e0,e^)  are  unit  vectors  in  the  spherical  system,  and 


m  e) 


1  —  ( 1  —  )  cos^  d 

y/l  —  (1  —  e‘^)  cos^ 


9{0;e) 


(1  —  )  sin  0  ce)s  6 

1  —  ( 1  —  e='' )  ce).s'‘^  9 


(30) 


(31) 


According  to  (16)-(17),  the  explicit  expression  fe)r  the  ce)e'fficieuits  A  9C  will 
be  of  the  form  of  finite  Legendre  transforms 


ft  =  {^+h)I  '  /I  TrtL 


9i  = 


Q  ^1-(1-€'‘)cos2  0 

r  (1  —  e^)sin0cos^ 


/ 

0 


i{£  +  1)  J  —  (1  —  cos^  9 


P^{cos9)fim^  9(19, 


(32) 

1/0=0,  (33) 


where  prime  indicates  differentiation  w.r.t.  (cos^). 

Both  integrals  vanish  for  odd  values  of  £.  It  is  shown  in  Appendix  B  that 
and  gi  can  be  represented  as  a  finite  sum  of  integrals  of  th('  tyi)o 


where  j  is  an  integer  and  cos^O  =  vl— This  integral  is  given  explicitly  by 
Gradshteyn  and  Ryzhik  (p.  131) 


^2; 


(2;  -  1)!!  f  TT 


Vj\ 


9  -  ^0  )  - 


sin  1^0 


COS^^  ^  i^O 


J-1 

7=1 


(2i-l)(2i-3)---(2i-2ry  +  l) 


2<I(j  -q) 


^0 


(35) 


Hence,  the  coefficients  /y  and  can  be  calculated  as  finite  doul^le  sums.  Once 
these  are  ready,  the  displacements  at  r  >  tq  are  a  special  ca.s('  of  {19)-(20). 

The  expressions  for  Up  and  Us  are  then  those  given  in  (21 )-  (22)  with  f  termi¬ 
nating  at  some  even  finite  value. 

Our  second  example  is  that  of  a  right  circular  cylindc;r  of  length  H  and  radius 
R.  Denoting  6q  —  tan~^  (^)’  normal  to  the  surface  is  given  by  n  =  ./e*,.  +geg, 
with 


j  COS 

9 

0  <9  <9ii 

/=  < 

sin 

9 

9o  <9  <  TT 

— 

—  cos 

V. 

9 

v  —  9q  <  9 

< 

TT 

—  sin 

9 

0  <  9  <  9o 

9=  < 

cos 

9 

9Q<9<n 

— 

[  sin 

9 

Tv  —  9q  <  9 

< 

TT 

(36) 


At  the  points  9  =  6q^  tt  —  6q  the  normal  is  not  defined.  According  to  (1G)-(17), 
the  coefficients  fi  and  gi  are  explicitly  given  by  the  integrals; 


fi  =  {2e  +  i) 


i+(-)' 


^0 

J  Pi{cos  6)  shiB  eoH  9(W 


7tf2 

+  J  Pi{cos  9)  sin^  9d9 1 , 


(37) 


do 


_  2£  +  i  ri  +  (-)^] 
i{e-\-i)[  2 


Oo 


J  P'^{cos  9)  9(19 

0 


7r/2 

—  J  P^{cos  9)  si-n?  9  cos  9 d9 
do 


0  cos  9 


(38) 
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Here,  again,  only  the  even  ^-values  contribute  to  the  sum.  It  is  also  shown  in 
Appendix  B  that  the  first  integrals  in  (37)  and(3S)  can  bo  evaluated  analytically, 
while  the  remaining  two  can  be  reduced  to  finite  suras  of  integrals  of  the  type  l2j  = 

cos^^  already  discussed  in  the  previous  example.  Therefore,  the  integrals 
of  and  Qi  reduce  here  too  to  an  evaluation  of  finite  double  sums  (Appendix  B). 


4.  The  Long-Wave  Approximation  and  the  Corrcspoiiding  Source 
Moment  Tensor 


In  cases  where  the  radius  of  the  cavity  is  smaller  than  tin.'  radiation’s  wave¬ 
length  such  that  kpTQ  <  1,  fcgro  <  1,  we  put 


1 


1 


-f 


1 


[x^+1  2(2t  - 


+  0 


(39) 


and  deduce  from  (8)  the  suitable  approximations  for  ('  =  0,2; 

■^0,1(2:)  =  -22x“^  -  - , 

■fb,2(a^)  =  -2x“^  +  ix“^  -1 - , 

F2,i  =  — 12zx~^  —  ix~^  , 

■^2,2(3^)  =  4Szx“^  -t-  3ix~^  +  •  ■  ■  . 


(40) 


■^0,3(2:)  =  2zx  ^  - 


la 


-^2,3(2;)  =  36ix  ^  i 


1  -  2cr 
1—5(7 


a-  ‘  -f 


1-2(7 


-b 


where  a  is  the  Poisson  ratio.  Using  these  results  wc  olitain  the'  approximations 


Ao  = 


4i 

x2y3 


A2  = 


12  7  —  5a 


x^i/^  1  -  a  ’ 

X  =  ksro,  y  =  kpTQ. 


(41) 


In  general,  for  ^  >  0 


~  -i 


{e  +  2) 


(2^-1)!!  m-3)\\ 


2x^+1 


+ 


(2(:-b  1)!! 


J-2 


+ 


(2('  +  l)!! 


2(r  - 1), 


(^  + l)(£  +  2) 

^1-2 


(2£-  1)!! 

x^+3 


+ 


1— i 

1 

Cd 

b 

1 

1 

(2f'-3)!!  ) 

b 

01 

1 

t— c 

(M 

i 

(42) 


+ 


(2£-M)!! 


^(£-1), 


60  - 


Ai{x,y)  ~  - 


(^  +  2)(2£-  l)!!(2£-3)!!  [(^'^  +  (  +  1)  -  a{2C  +  1) 


1  —  cr 


The  corresponding  values  of  and  [Eqs.  (25)-(2G)]  for  C  <  3,  arc 

A.  _  r  ■^0)2(^S^o)  _  i  r_  ru  . 


^0  =  /o- 


2A2  2A2 

=  -iU2  +  3g2)j^(^^yk„,-of,  (44) 

B2=[/2'  -52  Jj- 

=  i(f2  +  3g2)^^(^)a->nf'- 

When  this  approximation  is  applied  to  the  expression  of  the'  (jc.nc.val  {lisi)laceinent 
field  resulting  from  an  azimuthally- symmetric  surface  [Eqs.  (21)  (22)],  we  obtain 

^o.o( V) - (FTI)!! 


usinuj)  =  - 


f  ^[/Ke)  +  (£  +  iM^)]  1  /fcg- 
\(£2  +  £  +  1)_<t(2£  +  1) J  Vfcf. 

fpo(w)(l-(7)  -A  (fcsro)^"^^ 


L^^  Akpv), 


(2£  -  3)!! 


f=2,4,...  ^ 

fei^)  +  (^  +  l)g^(e) 

(£2  +  ^  +  l)_a(2£  +  l) 


N^Jhr). 


For  large  values  of  £,  the  dependence  of  /£  and  on  C  can  Ijc  estimated  from 
(16)-(17):  since  /  and  g  are  slowly  varying  functions  of  8  relative  to  P[{cos8)  and 
■^P^{cos  8)^  respectively,  it  can  be  taken  out  of  the  integrals.  Ic'ading  to 


TT 

e)  J  Pi{cos6)sin6(I0  =  2C  f\ 


~  /  -^Pi{cos  6)  sin  OdO  =  0  {C  evou). 


It  is  thus  guaranteed  that  for  ksVQ  <  1  contributions  from  suimnancls  in  (45)  for 
£  >  2  will  have  small  effects  on  both  Up  and  Us,  irrc.Hpcxtiva  of  tha  .specific  shape 
of  the  surface. 
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Note  that  if  instead  of  the  above  argument  wc  use  the  asymptotic  form  (Mag¬ 
nus  ei  al.^  1966) 


pp{cosO)  =  {-re'^-^/^- 


TT  sin0 


sin 


(47) 


(46). 


£>■1,  e  <  9  <  TV  — 
JL 


we  shall  obtain  fi  oc  (£  +  l)g'£  oc  with  a  faster  convergence  rate  than  in 


The  termination  of  the  sums  in  (45)  at  £  =  2,  yields  the  approximation 


fj.kp 

ipo(w) 


"■  “(/2  +  ^92) 


{k3rQf[f2  +  ^92]^— 2(^'.sO-). 


(48) 


fxks  '  '""'S-Ttr' 

Using  the  far-held  approximation  [(A-18)-(A-19)],  (48)  becomes 


{Ur}p  — 


fikp 


(fcpro)  <  -/o  +  3(./2  +  3(/2) 


Vpj 


7^(“=  ^-3 


ihjtV 


hr 


-  +3r/2)rr — ^  sin 261. 

Iflks  t  —  o<7 


(49) 

(50) 


It  is  useful  to  identify  the  displacement  field  in  (4S)  with  an  (ujuivalcmt  yoint  force 
system  operating  at  the  center  of  the  cavity  (Fig.  4).  This  can  easily  be  achieved 
when  we  write  down  the  displacement  fields  due  to  three  mutually  perpendicular 
dipoles  with  moments  {Mi  =  M2,  M3}  in  the  resp(!ctivc  :r,y,r  directions  (Ben- 
Menahem  and  Singh,  1981,  pp.  203-205): 


“3  =  [72(2X0,0  -  4X0,2)  -  4iVo,2] , 


U2  = 

Ui  = 


487rp 

iMikl 
48  TT/i 
iMik] 
487r/i 


[72(2X0,0  +  2X0,2  +  X2,2)  +  {'2Ni}^-2  +  -^2,2)]  ^ 
[72(2X0,0  +  2X0,2  -  X2,2)  +  2(iVo,2  -  -N2.2)] . 


(51) 


72  = 


2 


4 
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The  combined  field  {us  +  U2  +  }  can  be  split  into  the  following  coinpressional 

and  shear  fields 


Mb’-)  =  +  2Mi)^  L-„(l.;,r] 

o  47r/i  \Vp  J 


f- 


+  WV), 

ik^  ^ 

Us{ksr)  =  (Ml  -  M3) A^q,2(A:,s-?-)- 

A  complete  agreement  with  (48)  is  obtained  if  we  set  the  correspoudcnco 

/o  2(./2+3t/2)' 


3(1  -  a) 

Ml  =  M2  = - - - Wo 


I -2a 


I  —  OCT 


M 

Ms  =  - - - teo 


/o 


1  -  2a 


+  4 


(/2  +  3ty2) 


/  0(7 


(52) 


(53) 


4^  3  /  N 

^^0  =  y  i’oPo(w). 


Note  that 


<5M  =  Ml  -  Ms  =  -(/2  +  392)^^^-/^  «’o 

/  —  0(7 


>  0  line  source 
<  0  disc  source, 


(54) 


(55) 


cannot  assume  arbitrary  values  under  the  physical  (•t)iulitious  of  our  i)rol)lem.  To 
see  this,  let  us  examine  two  limiting  cases  of  special  interest: 


(I)  A  line-source  is  obtained  in  the  limit  k  =  I  for  the  spluuoidal  cavity  or 
^0=0  for  the  cylindrical  cavity.  We  assume  that  the  source  eiierg}'  (mo  is  fixed  such 
that  in  the  limiting  case  wq  =  limpo— 0°  [^^oPo]-  It  is  then  found  from  (32)-(33) 

Tq  —0 

and  (37)-(38)  that  for  both  the  spheroidal  and  cylindrical  cavities 


/0  =  J, 


h  =  -— 

32’ 


57r 
~  32’ 


Equatiems  (52)-(53)  then  yield  the  limiting  values 

21-k  297n  i 

■^3  =  =  0.23^0,  i\/i  =  M2  = ((>(),  for  (7  =  -, 

1 


Ms  =  0, 


A/i  =  M2  =^teo, 


for  a  =  -  . 
5 


(56) 


(57) 


(II)  A  disc-source  is  obtained  from  the  cylindrical  cavity  in  the  limit  6^0  =  f- 
Equations  (37)-(38)  then  render 


/o  = 


2’ 


<n  g . 


(58) 
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Ms  =  =  3.0Sii;o,  Mi  =  M2  =  =  O.luu’o,  for  cr  =  - , 

io4  io4  4 

M3  =  3ti;o,  Ml  =  M2  =  0  for  <7  =  -  . 

5 

Note  that  (55)  predicts  M2  =  Mi  —  M^,  under  the  condition 

/2  +  ^92  -  0.  (59) 

This  indeed  occurs  for  a  cylindrical  cavity  for  =  50.3°  {'2R  ^  1.2ff).  In  this 
case  the  cylinder  mimics  a  sphere,  with  no  radiated  shear  waves  (!),  independently 
of  the  medium’s  Poisson’s  ratio  [Figs.  5,  6]. 


Oblate  cavities 

As  long  as  the  observer  is  placed  in  an  unbouiuled  homogeneous  and  i.sotropic 
elastic  medium,  no  distinction  need  be  made  between  oblate*  and  prolate  cavity 
since  we  can  always  choose  the  2-axis  to  coincide  with  the  major  symnuitry  axis  of 
the  spheroid.  This  symmetry  is  broken  once  a  stress-fr('e  l)ouiidary  i.s  introduced, 
for  now  it  makes  a  great  difference  whether  the  cavity  is  aligned  normal  or  parallel 
to  the  free  boundary.  For  tunnels  with  axis  of  symmetrj’  i)arall(,'l  to  a  free  surface, 
the  radiation  patterns  can  be  obtained  by  rotating  the  ones  obtained  previously 
by  90°  relative  to  the  fixed  cartesian  coordinate  system  that  was  chosen  aljove.  In 
the  long-wave  approximation  this  amounts  to  an  eciuivalent  force  system  of  three 
orthogonal  dipoles  with  moments  {Mi,  M2  =  M3}.  Equation  (51 )  is  then  modified 
into 

“3  =  48ff  j*  -4£o,2)-4JVo,2], 

<J2  =  [72(2£o,0  +  2Lo,2  +  £2,2)  +  (2A'„,2  +  )] ,  (60) 

zAjTi  Jc^  ^ 

[72(2X0,0  +  2Lo,2  -  ^2,2)  +  2(.No,2  -  -^^2,2)]  • 

In  this  new  configuration  the  x-axis  is  aligned  along  the*  smaller  moment  Mj,  now 
parallel  to  the  free  surface.  The  combined  field  {t7;j  -|-  U2  -f-  a  I }  is  thc'n  expressible 
as 

Up{kpr) 


Us(kar) 


+  £„-,o(V) 


+  (M3  -  Ml) 


ikg  /  Vs 


(61) 


^2,2  (/■■;/'•  )^ 


247r/i  \vp 
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The  source  moment  tensor 

Every  c^olax  source  can  be  represented  by  a  second  order  symmetric  carte¬ 
sian  tensor  M  (Ben-Menahem  and  Singh,  1981,  pp.  lGS-171).  Its  corresponding 
displacement  field  is  given  by  the  expression 

u(r)  =  M;grado{G(r  1  (61a) 

i f 

where  G  is  the  Green’s  tensor  of  the  elcistic  medium.  Wo  ha.vo  first  shown  that  the 
displacement  field  of  a  symmetrical  cavity  with  a  syiiimctiy  axis  in  the  ^-direction 
is  equivalent  to  the  combined  fields  of  three  dipole.^,  each  along  a  coordinate  cixis, 
at  the  center  of  the  cavity,  with  the  respective  strengths  M 1 1  =  i\/| ,  M22  ~  -^^2  = 
-^1?  =  M3. 

Thus,  in  the  light  of  (53),  the  moment  tensor  of  tlu'  cavity  can  be  written  in 
the  dyadic  form 

M  =  Ml  I  +  (M3  —  A'l]  )c;jC3  (Clt) 

< '  ^ 

where  M,  I  is  known  as  the  isotropic  part  of  the  tensor,  and  —  M]  )c3e3  is 
the  deviatoric  part.  This  last  part  is  the  true  signature  of  the  iioii-sphcrical  cavity 
which  makes  it  distinct  from  both  symmetrical  explosions  and  (.'arthcpiakes. 


5.  Surface  Waves  from  Explosions  in  An  UndcrgrouiKl  Cavity 


Knowing  the  source-fields  in  an  unbounded  domain,  the  surface-wave  fields 
in  the  presence  of  a  stress-free  planar  boundarj'  can  lx*  (U'aluatc'd  in  a  routine 
manner.  This  becomes  necessary  if  one  wishes  to  calculate  theoretical  waveforms 
of  Rayleigh  waves  produced  by  explosions  in  underground  tunnc'ls.  Since  usually, 
the  long  axis  of  the  cavity  is  parallel  to  the  free  surface,  we  must  use  as  our  source- 
field  the  displacement  produced  by  an  oblate  cavity  as  given  in  E(j.  (61),  where 
we  now  denote  it  as 

=  Up{kpr)  -f-  Us(/.:s/  ).  (62) 

It  is  convenient  to  express  the  l.h.s.  of  (61)  in  terms  of  the  cyliudrical  vector 
harmonics  via  (A-l)-(A-8).  It  yields  the  source  disi)lac('m('nts  and  associated 
tractions  in  the  form 


where 


oo 


40) 
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(/.•) 


40) 


,  40)^- 

1  -  tn  /] 


(63) 


(64) 
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(^,m)  =  (0,0),  (0,2),  (2,2] 


(g)  I-I 


a-2  -  A:2 


'k^  -  k'i 


T]s  =  Z€- 


e  =  Sgii(c  -  '()),  (66) 


with  the  coefficients  and  defined  in  (25)-(20). 

The  stress- vector  (normal  traction)  associated  with  the  source  displacements, 
is 

e,  ■  =  e^Adivu^®)  -t-  +  (7(‘’)V). 


=po(u^)  Jyi 

0 


m-C+lfW  I 


with  the  dimensionless  partial  stress  expanded  as: 

P  _i_  p 

Mm  —  rn  'T  i  m  “r  ^  tii  ^ 


LyCO)  _  a^(0)  , 

Ki  m  ^  AJXm  “T 

7,7(0)  _  dz^m^ 


/.rrtO) 

IJin 


Using  (65),  this  becomes 
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A'S’  =  2A,pp(n,) 


2fc2  _  kl 


-4Be(^ 


=  -4eA(P, "•{’Ip  )«-■'' 


Z/p/u 

P/ 

—  I'vlz 


idijs 


+  2B,|i 


P  /'o7,2  7,2’, 

437(2^ 


2i!>  =  -2eB<|£ 

tZp 


p  L 

l2 


c/pr(v.)  - 


idi}s 


(70) 


(£,m)  =  (0,0),  (0,2),  (2,2). 

^te  that  in  our  formalism,  }  are  diincnsioiilc'ss,  but 

T(u^  )  j  have  the  respective  dimensions  of  displacement  and  stress. 

We  can  now  state  our  boundary  value  problem  in  terms  of  tlu'se  source  dis¬ 
placements  and  stresses  and  the  yet  undetermined  corrc^spondiiifi;;  entities  of  the 
half-space.  Assume  an  expansion  of  latter  in  terms  of  cylindrical  vc'c’tor  harmonics 
[Appendix  A], 

Um  =  ^mPm  +  ymPrn  +  ^u^C^n- 
Tffi  =  A  772  -\-Y Bfxi  +  Z  C\]i , 

where 

^  “b  ^ 

Vm  —  drn^^  ^  d*  ^  ? 

Here,  (^^m^  C772)  ^-re  undetermined  amplitude  coefficients  (^f  dimensions 

{length  Jength,  1)  respectively 

kXrn  =  {2k^  -  kl)ame^^^  +  2kub,nc‘'^^\ 

kYm  =  2upkame^r>^  +  {2k^  -  (73) 

kZfn  ~  Cm C  ^  . 

Again,  {^m,  yrn?  A777,  V  772,^777}  are  dimensionless  pln'sical  entities.  The  total 
field  is  u  +  Applying  the  condition  of  a  free  surface  at  r  =  h 

T[u  +  u(^)]  =  0  at  .:  =  /?,  (74) 


(71) 


(72) 


67 


we  get 


=  yS^^  +  rm  =  o,  +  =  o  at  -  =  /?.  (75) 

Solving  for  the  coefficients  {aTn,bmiCm]  and  substituting  the  results  in  the  ex¬ 
pression  for  the  total  displacement  field,  we  find  for  the  field  on  the  free  surface 
(e  =  1,  z  =  h) 


^total  ^  -(0)  ^ 


Po(w) 


/£ . 

0  ^ 


u  =  — - — -  /  >  i"~  ''  '  * u, 

fikp 


f-fl  -lolal^^ 
9/.2  ’ 
-"■p 


=  UpPm  +  UeBm  +  C/cCnt, 


(76) 

(77) 


Up  =  +  2B,h 

n>rt 


_ i 

■P  L^s^s  J 


idiis 


04  pm/  -i/p/i  op  dPl'^Tjs)  P  Pj,  Rtik)  , 

--2AtPt  W-Rjty'  ' 


Us  =  4^,(  f  )Pr(,/)e-^''‘  -  2B,|  (^) 


2^(-P(  W  TFi?(4)"  ' 


(^,m)  =  (0,0),  (0,2),  (2,2), 

(78) 

R{k)  =  (2k^  -  k^f  -  AJPupUs. 

Rp{k)  =  {2k‘^  -  k'lf  -b  AlPupUs. 

(79) 

i7l(^)  =  R-^{k)  -  ^k^{2kP  -  k^),  R2{k)  =  R+{k)  -  AupUsi'llP  -  k^). 


The  source-terms  in  the  expressions  for  Uq,  Up  and  Up  do  not  coiitril^ute  to  the 
surface- wave  field  and  can  therefore  be  ignored  in  this  study.  The  remaining  terms, 
which  have  R{k)  in  their  denominator  [Eqs.  (62)-(G3)],  contribute  to  the  Rayleigh 
wave  through  the  pole  of  R{k)  =  0  in  the  lower  complex  /.--plane  [Re  w  >  0;  see 
e.g.  Ben-Menahem  and  Singh,  pp.  263-265]. 

The  relevant  integral  in  (60)  has  the  form  ./„,(/<: A) 

=  ^  Um\kA)^^dk,  where  Hm\kA)  are  the  Hankel  hmetions  of  order  m 
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of  the  second  kind.  The  pole  of  R{k)  =  0  is  at  Icji  =  7A\^.  wIiotc'  7  is  a  root  of 
the  equation  (27^  -  l)^  =  4^2 _  i  ^.^2  _  ^  ^  ^  j  ^  = 

(2-  =  1.0875. 

The  residue  of  the  above  integral  a.t  k  =  kji  is  given  by  {  Hl,i\kjiA) 

^  =  2.51220  .... 

The  terms  in  (62)-(64)  which  do  not  have  R{k  )  in  their  dciioniinator  (direct 
waves),  do  not  contribute  to  the  residue  field  and  hence  do  not  g(;n(;rate  surface 
waves.  Inserting  in  (76)  all  the  relevant  entities  from  Eqs.  (7S)-(79),  using  the 
relations  kji  =  jks,  ^  (^  =  ?)>  7p  =  -  1,  'h  =  -  1  nnd 


Rlikji)  =  -2kti2-f^  -  1),  R2{k)  =  l4('^ - 


Hii\kjiA) 


e  2  ^  -f- 0(A7^ A) 


irkjiA 


the  components  of  the  spectral  Rayleigh-wavo  displac('iu('iits  for  C  <  4  at 
(z  =  /i,  A,(/?)  which  fall  off  with  distance  like  arc 


5i  —  —  [(2Ao  +  A2)  —  97^ A2  cos*^  , 

^2  =  -“(27^  —  1)52  cos^  if. 


6.  Numerical  Results  for  Radiation  Patterns  of  Body  and  Surface 
Waves 


For  the  numerical  calculations  of  radiation  patterns  of  l)0(ly  waves  generated 
by  explosions  in  prolate  spheroidal  and  cylindrical  ca.viti('s,  we  used  Eqs.  (21)  and 
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(22)  normalized  by  We  calculated  the  first  five  terms  in  the  infinite  series 

(i.e.  expansion  up  to  the  8^^  order  in  spherical  harmonics).  All  the  vector  spherical 
harmonics  Lqi  and  Nq(^  (€  =  0, . . . ,  8)  were  calculated  exactly  using  the  formulas 
given  in  Appendix  A.  The  functions  x,  ^^,2)  -ft,;]  calculated  using 

formulas  (8).  The  shape-dependent  coefficients  f[  and  (j[  were  calculated  exactly 
using  the  equations  (32)-(33)  for  prolate  spheroidal  cavities  and  (37)-(3S)  for 
cylindrical  cavities.  The  integrals  in  Eqs.  (32)-(33)  and  (3T)-(3S)  were  calculated 
as  finite  double  sums,  using  the  analysis  given  in  Appendix  B. 

Spheroidal  cavity 

Figure  5  presents  radiation  patterns  of  body  wat’cs  in  the  c-A  plane  for  explo¬ 
sions  in  prolate  spheroidal  cavities  (the  z-axis  corresponds  to  the  axis  of  sjunmetry 
of  the  cavity).  Radiation  patterns  were  calculated  for  =  0.01  (long  waves), 
kpT  =  100.0  (far-held).  The  six  parts  of  Figure  5  correspond  to  the  values  of  the 
aspect  ratio  e  =  ^  =  1.0  (a  sphere),  0.8,  0.668,  0.414,  0.199  and  0.0  (line  source). 
For  a  spherical  cavity  (Fig.  51)  the  radiation  pattern  of  P-waA‘c.s  is  symmetric  and 
no  S-waves  are  generated.  Since  in  the  far-held  terms  of  orcU-r  are  lu^gligible, 

the  P-patterns  are  those  of  the  radial  component  of  th(’  disi)lacement  vector  (23), 
while  the  S-patterns  are  those  of  collatitudinal  comx^onent  of  the  displacement 
vector  (24). 

As  the  shape  of  the  cavity  deviates  from  sphcnical,  tin.'  P-wave  radiation  pat¬ 
tern  becomes  stretched  in  the  direction  perpendicular  to  the  major  (the  longest) 
axis  of  the  cavity.  At  the  same  time  shear  waves  arc.'  geneiaitc-d. 

One  can  see  that  radiation  patterns  are  dominated  by  Loo,  -^02  foi’  P-waves 
and  Nq2  for  S-waves  [the  calculated  coefficients  of  highc'r  (C  >  2)  order  spherical 
harmonics  were  several  orders  of  magnitude  smaller  for  k],i\)  1].  This  is  consis¬ 

tent  with  the  conclusions  of  Section  4,  that  the  long  wave  radiation  patterns  are 
dominated  by  dipole  terms.  As  the  shape  of  the  cavity  is  stretched  from  a  sphere 
to  a  line,  the  vertical  (along  z-axis)  dipole  becomes  smaller,  and  the  horizontal 
(along  X-  and  y-axes)  dipoles  become  bigger.  Therefore  the  P-wavc;  radiation  pat¬ 
tern  becomes  stretched  in  A-direction,  and  shear*  waves  are  generated. 

Cylindrical  cavity 

Figure  6  presents  radiation  patterns  of  body  wavc-s  in  the*  z  A  plane'  for  a 
cylindrical  cavity  (the  z-axis  corresponds  to  the  axis  of  symmetry  of  the  cavity). 
Radiation  patterns  were  calculated  for  kpVQ  =  0.01  (long  wavers).  kpV  =  100.0 
(fax-field).  The  six  parts  of  Fig.  6  correspond  to  the  cylinder  ai)pert\ire  angles  of 
^0  =  0.0  (line  source),  0.8779  (50.3°),  ^  and  ^  (elisc  source).  As  for  the 

previous  case  of  the  prolate  spherical  cavity,  the  radiation  x)atrerns  are  dominated 
by  dipole  terms. 
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For  ^0  =  0.8779  (50.3°)  all  three  dipoles  Mi.  d/2,  d/.-j  have  equal  values  [i.e. 
fl  +  ^g2  =  0  in  equations  (53)]. 

In  order  to  study  the  effects  of  higher  order  (/  >  2)  .spherical  harmonics, 
we  calculated  radiation  patterns  for  wavelengths  comparable  with  the  size  of  the 
cavity  (the  radius  of  the  equivalent  sphere).  Figure  7  presents  radiation  patterns 
of  body  waves  calculated  for  kpVQ  =  1.0,  kpV  =  10*.  One  can  see  that  higher 
harmonics  change  only  a  little  the  radiation  patterns  for  this  wavelength  [probably 
the  strongest  effect  can  be  seen  in  Fig.  6IV  for  6^0  =  0.S779  (50.3°).  The  S-wave 
radiation  pattern  is  given  primarily  by  Nq^]. 

Higher  order  (/  >  2)  spherical  harmonics  play  a  significant  rok;  only  for  wave¬ 
lengths  much  smaller  than  the  size  of  the  cavity  {kpvo  >1).  However,  for  realistic 
cavity  sizes  (100  meters  at  most)  this  wavelength  are  not  observed  in  far-field 
seismograms  due  to  attenuation.  Therefore  we  may  conclude  that  the  far-field 
radiation  patterns  of  body  waves  generated  by  explosions  in  arbitrary  cavities  will 
be  dominated  by  dipole  terms  (i.e.  two-lobe  radiation  pattern  for  P-waves  and 
4-lobe  patterns  for  S-waves). 

This  conclusion  is  consistent  with  the  results  of  dirc'ct,  numerical  simulations 
by  Stevens  et  al.  (1991)  who  presented  radiation  pattern  for  an  explosion  in  a 
prolate  spheroidal  cavity  with  aspect  ratio  e  =  0.25. 

Figure  8  exhibits  the  spectral  amplitude  depend(uic(;  of  tlu'  field  on  the  apper- 
ture  of  the  cylindrical  cavity.  We  note  a  node  for  S-wavc's  at  B{)  %  50.3°,  at  which 
the  P-waves  pattern  is  spherical.  The  S/P  energy  ratio  is  extremal  at  Bq  =  0  (line 
source)  and  ^0  =  f  (disc  source),  but  is  small  ovt'r  the  wide  plateau  (9q  =  | 

Corner- frequencies 

In  order  to  calculate  the  frequency  dependence  of  body  wave's  amplitudes  we 
had  to  make  an  assumption  about  the  time  dependence  of  the  pr(!ssurc  pulse  P{t) 
applied  at  the  surface  of  the  cavity. 

The  most  straightforward  assumption  is  that  P{t)  is  given  l;)y  the  step  function 


P{t)  = 


1  t>0 

0  t  <0. 


(81) 


the  Fourier  transform  of  which  is  given  by 

1 


■Po(^)  =  —  ^  0). 

tU) 


(82) 


We  calculated  amplitudes  of  body  waves  given  bv  (21)  and  (22)  using  (82)  for 
^  =  10"^. 
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Figure  9  presents  frequency  dependence  of  P-\vave  amplitude  generated  by  the 
pressure  pulse  (81)  applied  at  the  surface  of  cylindrical  ca,vity.  Figure  9  corresponds 
to  cylinders  with  apperture  angles  =  0.0  (line  source),  q,  0.8779  (50.3°),  ^ 
and  Y  (disc  source). 

Each  of  the  figures  contains  amplitude-frequency  dependences  for  signals  ob¬ 
served  at  the  5  collatitude  angles  ^  =  0,  ^  and  y-  We  plotted  the  logarithm 

of  the  amplitude  versus  the  logarithm  of  the  dimensionless  frequency: 

t>ro  =  (83) 

Vp 


All  calculations  are  made  for  ^  =  10 

r 

It  is  clearly  seen  that  the  “corner-frequency”  in  all  these  plots  corresponds  to 
the  wavelength  equal  to  the  radius  of  the  equivalent  sphere  of  the  cavity. 

Figure  10  presents  the  amplitude-frequency  dcpeiuU'ncc's  of  S-wa.ves  for  the 
same  situation.  The  plots  contain  amplitude-frecpu'iicy  dc'pc'iidencc's  of  signals 
observed  at  collatitude  angles  ^  =  f  >  f  and  ^  (the  ainplitude  of  S-waves  radiated 
in  the  direction  along  and  perpendicular  to  the  axis  of  s3unmetry  of  the  cavity  is 
negligible  in  comparison  to  the  P-wave  amplitude). 

For  these  plots  we  used  dimensional  frequency'  given  bj" 


ksTQ  = 


rpu; 


(84) 


It  is  again  obvious  that  the  “corner-frequency”  corresponds  to  a  wavedength  equal 
to  the  radius  of  the  equivalent  sphere. 

We  may  therefore  conclude  that  the  size  of  the  cavitj'  wherein  the  explosion 
takes  place  can  be  deduced  from  the  corner-frequ(>iicy,  •providc.d  that  there  was  full 
decoupling  and  that  the  shock  wave  —  medium  interaction  occurred  in  tlic  elastic 
regime. 


Surface  waves 

For  numerical  calculations  of  surface-wave  radiation  paf.tc'rns  we  used  the 
dipole  approximation  derived  in  Section  5. 

Figure  11  presents  radiation  patterns  of  Rajdeigh  surface  wave's  geiK'rated  by 
explosions  in  cylindrical  cavities  with  an  horizontal  axis  of  sj'inmetrj’  parallel  to 
the  free  surface.  Radiation  patterns  were  calculatc'd  for  hpi'o  =  0.01  (long  wave), 
kpr  =  100.0  (far-field),  kpr  =  0.0  (shallow  source),  i.c'.  for  a  s(.)urce  at  depth  much 
smaller  than  the  wavelength.  The  pattern  in  Fig.  11  correspond  to  cylindrical 
cavities  with  apperture  angles  =  0.0  (line  source),  f ,  f ,  0.8779  (50.3°).  ^  and 
^  (disc  source). 
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One  notices  that  for  a  source-depth  much  suialler  than  the  wavelength,  the 
Rayleigh- wave  radiation  patterns  resemble  the  ones  for  the  corresponding  P-waves. 

Figure  12  presents  radiation  patterns  for  the  same  explosions  ljut  for  a.  "deep 
source  ’  with  kpr  =  1.0.  At  this  depth,  the  radiation  patterns  change  considcu-ably 
relative  to  the  corresponding  ones  at  h  =  0,  due  to  the  amplitude  decay  with 
depth. 

Corner  frequencies  for  surface  waves 

Figure  13  presents  amplitude-frequency  dependences  of  Raykagh  surface  waves 
generated  by  pressure  pulse  (81)  applied  at  the  surface  of  cylindrical  cavities  with 
horizontal  axes  of  symmetry.  The  amplitude-frccpiency  dependences  tue  calcu¬ 
lated  for  kph  =  0.0  (shallow  source),  ^  =  10“'^  and  plotted  against  dimensionless 
frequency. 

The  various  patterns  in  Fig.  13  correspond  to  caviti(;s  with  appenture  angles 
^0  =  0.0  (line  source),  f ,  f ,  0.8779  (50.3°),  ^  and  f  (disc  source).  The  five 
curves  at  each  figure  correspond  to  the  singal  observed  at  tlie  5  azimuthal  angles 
^^  =  0’ 

The  left  side  of  all  the  amplitude-frequency  dopcndciiccs  ha.ve  a  slope  of  ^ 
(compare  with  the  same  curves  for  body  waves  for  which  the  l(‘ft  hand  side  is 
horizontal).  This  is  because  the  surface  waves  decay  at  - — --777  ratlier  than  7-^ 

(body  waves),  which  brings  in  an  extra  factor  of 

For  a  “shallow  source”  the  “corner- frequency”  corresponds  to  a.  wavelength 
that  is  approximately  equal  to  the  radius  of  the  equivalent  splu'rc'. 

Figure  14  presents  the  same  amplitude-frequency  dependences  but  for  a  "deep 
source  with  kph  =  1,0.  One  can  see  that  the  ainplitude-frecjiuncy  dependences 
are  mainly  controlled  by  the  exponential  terms  and 


7.  Discussion 

The  radiation  of  seismic  body  and  surface  waves  from  ('xplosions  in  under¬ 
ground  tunnels  has  been  calculated  in  the  far-field  of  tlu'  ('la.sti(‘  zoiu'.  The  following 
assumptions  and  approximations  were  made: 

•  The  reflected  body-wave  field  from  the  free  surfa.cc'  was  ignored;  onh^  direct 
P  and  S  signals  were  considered.  Since  these  wave's  are'  trulj"  diagnostic  of 
the  source,  it  is  perhaps  preferred  to  isolate  the  direct  fic'ld  from  the  data 
and  compare  the  Fourier  transformed  P  and  S  wavc'forms  with  our  theoretical 
model. 
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Conversion  of  body-waves  to  Rayleigh  waves  and  vice  versa  was  ignored. 

Our  theory  is  valid  only  for  wavelengths  which  are  much  larger  than  the 
dimensions  of  the  cavity.  Shorter  wavelengths  are  assumed  to  he  aljsorbed 
and  not  present  in  the  data. 


It  is  clear  from  Figs.  5,  6  and  7  that  the  radiation  pattern  for  P-wa.vcs  is  always 
dipolar,  i.e.  it  can  be  represented  by  a  single  equivalent  dipole  with  at  most  two 
lobes  and  one  nodal  line.  The  S-wave,  on  the  other  hand,  is  quadTupolar^  winch  is 
the  resultend  of  two  perpendicular  dipoles. 


We  know,  however,  that  most  earthquakes  have  cpiaxlrupolar  P-wave  patterns 
(Frohiich,  1994)  which  arise  from  the  double-couple  nature  of  its  (‘ciuivalcnt  force- 
system.  Thus,  for  the  sake  of  comparison  with  (Gla),  the  monicnt  tensor  of  a 
standard  earthquake  hcis  the  dyadic  form 


M  =  —  6262) 

=  Mo  I  —  Mo(2e2^2  +  ). 


(S5) 


By  removing  the  isotropic  part  from  (Gla)  and  (85),  oiu*  can  (^xainiiu'  the  deviatoric 
part  of  the  cavity  explosion  to  see  whether  it  conforms  to  a.  standard  (‘arthquake 
source.  Obviously  the  two  deriatoric  parts  are  worlds  appart  l)('caus('  of  the  extra 
term  {— 2Moe2^2}*  Thus,  the  explosion  will  look  very  strongly  like  a.  ^-non  double 
couple  earthquakeh  This  can,  and  should,  serve  as  a  ck'finite  t('st  for  a  suspicious 
clandestine  explosion. 


Stump  ei  aL  (1994)  analyzed  the  near-source  sesimograms  recorded  from  the 
Coalora  nuclear  explosion  (Yucca  Flats,  Nevada  Test  Site,  February  11,  1983). 
They  showed  that  the  source  moment  tensor  had  an  isotropic  part  5  to  10  times 
greater  than  the  deviatoric  part,  which  according  to  thenr  analysis  was  mostly  due 
to  spall.  We  suggest  that  the  data  inversion  scheme  used  l>y  the  al)ove  authors  can 
be  applied  to  signals  arising  from  underground  explosions  in  order  to  determine 
the  influence  of  the  shape  of  the  cavity  upon  the  deviatoric  component  of  the 
source’s  moment  tensor. 


The  theory  is  valid  for  a  homogeneous  elastic  nu'dium.  However,  known 
algorithms  for  multilayered  media  can  be  applied  to  our  source  fields,  to  generalize 
the  results  for  multilayered  media  that  takes  int(j  account  such  phenomena  as 
dispersion,  attenuation  and  scattering.  The  fast  convergence  of  our  computational 
scheme  will  not  be  affected  by  this  generalization,  since'  this  conve'rge'iice  depends 
only  on  the  shape  of  the  cavity  and  not  on  the  propc.i-tie's  of  rlu'  elastic  medium. 

On  the  other  hand,  the  propagation  of  the  source  fields  through  a  vertically 
heterogeneous  media  will  generate  SH  and  Love  waves,  the  analysis  of  which  may 
render  additional  information  on  the  shape  of  the  cavity  and  the  spc'ctrum  of  the 
initial  pulse.  Indeed,  Love  waves  and  Lg  signals  wc're  oljscrved  from  underground 
decoupled  nuclear  explosions. 
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The  present  paper  does  not  address  directly  the  problem  of  seismic  deco  upling, 
namely,  the  dependence  of  the  signals  spectra  upon  the  shape  and  strength  of  the 
initial  pressure  pulse,  and  the  dependence  of  the  seismic  yield  on  the  shape  of 
the  cavity  and  the  rock-mechanical  properties  in  the  neighborhood  of  the  cavities 
boundary. 

The  main  thrust  of  the  present  paper  was  concentrated  on  the  extension  of 
Sharpe’s  naive  model  to  non-spherical  cavities.  The  suljject  of  decoupling  and  its 
implication  to  monitoring  of  small  nuclear  explosions  in  underground  cavities  will 
be  treated  in  a  sequel  paper  in  this  Journal. 
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Legend  of  Figures 


The  geometry  of  explosion  in  cavities:  (1)  spherical  cavity  in  an  un¬ 
bounded  solid;  (II)  prolate  spheroidal  cavity  and  its  •ccpiivalent  si)here’; 
(III)  oblate  spheroidal  cavity  in  a  half-space  with  a  free  surface  =  h). 


Major  axis  cross-section  of  a  prolate  spheroidal  cavity  r{0) 
2 

^  =  1  —  (f )  .  The  equi- volume  condition  renders  tq  =  a 


\/\—k  COS^  d  ’ 

for  the 


equivalent  sphere. 


Cylindrical  cavity  of  radius  R  and  height  H.  The  parameter  Oq  = 
^9  controls  the  radiation-pattern  of  the  source. 

Simulation  of  cavities  by  a  combination  of  three  mutually  orthogonal 
dipoles:  (a)  vertical  cavity  with  Mi  =  M2  >  (b)  horizontal  cavity 

with  M3  =  M2  >  Ml- 


Radiation  patterns  of  body  waves  displaceinents  gc'uc'ra  tc'd  l)y  explosions 
in  prolate  spheroidal  cavities  with  different  asi)ect  ratio  f  =  y  (Fig.  2). 
Patterns  are  drawn  in  a  vertical  symmetry  plane  embedding  the  r-axis, 
with  the  collatitude  angle  6  increasing  from  zero  {z  =  0)  clockwise:  (I)  e  = 
1.0  (sphere);  (II)  e  =  0.8;  (III)  e  =  0.668;  (IV)  e  =  0.414;  (V)  e  =  0.199; 
(VI)  e  =  0.0  (line  source).  Solid  line  —  P-wa.ves  (radial,  a,.);  dashed  line 
—  S-waves  (collatitudinal,  ug).  kpro  =  0.01  (long  wave);  =  100.0 
(far-held);  cr  = 

Radiation  patterns  of  body  waves  generated  by  explosi(jns  in  cylindrical 
cavities  with  different  apperture  angles  (I)  ^0  =  0.0  (line  source); 
(II)  00  =  I ;  (III)^  00  =  f ;  (IV)  00  =  0.8779;  (V)  0o  =  (VI)  =  f 
(disc  source).  Solid  line  —  P-waves;  dashed  line  —  S-wav('s.  kpro  =  0.01 
(long  wave);  kpr  =  100.0  (far-held);  cr  =  ^. 

Radiation  patterns  of  body  waves  generated  Ijy  explosions  in  cylindri¬ 
cal  cavities  with  different  apperture  angles:  (I)  0o  =  0.0  (line  source); 
,.(11)  ^0  =  f ;  (HI)  ^0  =  I ;  (IV)  ^0  =  0.8779;  (V)  ^0  =  (^L  ^0  =  f 

(disc  source).  Solid  line  —  P-waves;  dashed  line  S-waves.  KpVo  =  1.0; 
kpr  =  10^;  <T  =  5. 

Dependence  of  the  body  wave  amplitude  on  the  apperture  angle  for  a 
cylindrical  cavity.  Observation  angles:  - 0°;  “  f ! - f; . 

X; - f-  (I)  |Fp|  versus  20o/r-,  (II)  |(7.s|  vctsus  20o/n\  (III) 

versus  26o/7r,  Note  the  nodal  point  (N)  at  which  the  P-ia.diation  patterns 
are  spherical  (no  S-motion). 


Amplitude  (ln|up|)  frequency  (ln[/:pro])  dc'peudc'ucx's  of  P-waves  gener- 
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ated  by  a  pressure-step  applied  at  the  surface  of  a  finite  cylindrical  cavity 

observed  at  different  collatitude  angles:  - =  0:  . 

^  =  f; . ^  =  - ^  =  j.  Apperture  angles  of  cylindrical  cavi¬ 
ties  assume  the  values:  (I)  =  0  (line  source);  (II)  6^0  =  |;  (III)  =  j; 

(IV)  00  =  0.8779;  (V)  0o  =  x'-  (^1)  ^0  =  f  (disc  source).  ^  =  10"^ 
cr  =  ^.  Note  the  corner-frequency  at  kpVQ  =  1. 

Fig.  10  Same  as  8,  but  for  S-waves.  Observation  angles  are:  -  9  =  - 

Q  _  TT  .  a  _  Stt 

. ^  =  -g-' 

Fig.  11  Azimuthal  radiation  patterns  of  the  vertical  Rayleigh  surface  wave  dis¬ 
placement  generated  by  explosions  in  cylindrical  cavities  with  different 
apperture  angles:  (I)  0o  =  0  (line  source);  (II)  0o  =  f ;  (HI)  0o  =  f : 
(IV)  00  =  0.8779;  (V)  0o  =  ^;  (VI)  00  =  ^  (disc  sourc(').  The  r-axis  of 
the  cavity  (Fig.  3)  is  parallel  to  the  free  stirfacc'.  Note  that  the  pattern 
of  body  waves  and  surface  waves  from  the  horizontally  cylindrical  ca.vity 
are  both  in  the  same  horizontal  plane ^  and  therefore  c:in  bt'  compared. 

Fig.  12  Same  as  10  but  for  a  deep  source  {kph  =  1.0). 

Fig.  13  Amplitude  (In  luzl)  frequency  (ln[fcpro])  dependences  of  Rtiyleigh  surface 
waves  generated  by  a  pressure  step  applied  :it  the  surfivee  of  cylindrical 

cavity,  observed  at  different  azimuthal  angles:  - ■k  =  0; - F  = 

- ip  =  . p  z=  - p  =  I.  Apperture  tuigles  of  cylindrical 

cavities  are:  (I)  0o  =  0  (line  source);  (II)  0o  =  |;  (III)  0q  =  ^;  (IV)  0o  = 
0.8779;  (V)  0o  =  ^;  (VI)  0o  =  ^  (disc  .source).  ^  =  10“'*;  cr  =  |; 
kph  =  0.0  (shallow  source).  Note  the  corncT-fraquency  at  /.p/'o  =  1. 

Fig.  14  Same  as  12  but  for  deep  source  kph  =  1.0. 
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Appendix  A:  Fundamental  Elastodynamic  Vectors.  Associated 
Functions  and  Coordinate  Transformation  Relations 


The  Hansen  elastodynamic  eigenvectors  in  spherical  ccKnclinatcs  (r,6,y)  with 
unit  vectors  (er,e^,e<^)  are: 


^ml  ~  +  V +  1)  ^mt 


-  zf(kpr) 


kpr  ’ 


^ztik.r) 


+  VI(e  +  T)Bmi 


k.r 


(-4-1) 


ks  = 


=  V^{^  +  l)C^izf{ksr). 

The  spherical  vector  harmonics  are: 

Pmi  =  erPY^{cose)e^^^, 


•jj 


Cv 


1  d 


ppicossy 


(-4-2) 


iriiKp 


In  (A-l)-(A-2),  z^(x)  =  ji{x),  spherical  Bessel  functiou  of  the  first  kind  and 

z^  (x)  =  Rx),  spherical  Bessel  function  of  the  second  kind.  P"'(cost/)  is  the 
associated  Legendre  polynomial  of  degree  ^  and  order  in.  A  prime  {')  denotes 
differentiation  w.r.t.  the  argument,  and  ki  =  kp  or  ks. 

In  circular  cylindrical  coordinates  (A,  (^,2)  with  unit- vectors  {c^.tVVs},  the 
Hansen  vectors  assume  the  form: 


Lt  =  ■^{±UpPm  +  kBm)e^‘'r  =  , 


Mm  = 


(-4-3) 


where 


Up  =  ^Jk^-kl  Us  =  ^/k^-kl  Ym{kA,  VV  =  .7,„  (/.^A)r 


PmikA.if)  =  e^YmikA.ip),  m  =  0,1,2 . 


Bm{kA,y?)  -  (  ^A-g^  + 


1  d 


Cm{kA,(p)  —  (  e^-:—r— —  —  e 


kA  d<p 
^  d 

kA  dip  '^'^Wa 


Yui[kA.p>), 

):  ,n(kA.  pi). 


(-4-4) 


-  81' 


We  shall  need  the  results, 


div  Pjn  =  0,  Bz  div  Bm  =  -kPm,  div  =  0, 


{A -5) 


■  (Vi^ rzi  "H  ■PmV)  —  kBiti, 

ez-{VBm  +  Bmy)  =  0,  (*4-6) 

B-z  •  (\^Cm  “h  C'mV)  =  0. 

Transformation  of  spherical  to  cylindrical  eigenvectors  was  treated  Ijy  Bcn-Mcna- 
hem  and  Singh  (1981,  p.  78).  It  is  shown  there  that  the  exterior  Haiiseu  vectors 
in  spherical  and  cylindrical  coordinates  are  linked  through  the  integral  relations. 


=  —  * 


_  1  -m-i+l 


C50 

|rs;>(ipA)pr(,p) 


kdk 


I  N!n\k.A) 

^  0 
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oo 


\dpri’h} 


l-^dk 
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m—^+1 


dris  J  -  I'i 

7  — (Ik 

J  iksA)PrOu)-=f=  . 


(.4-7) 


OO 


dproh}] 


khlk 


oo 


s(0) 


_,m-e+l  j  (k.A)P^(r,, 
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^Vs  J  \/}x.^  — 

dk 


Here 


e  =  sgn(2  -  20 ),  Vp  =  ze 


_  • 


Vs  =  <(- 


h 


N^\ksA)  —  —(kPm  ~  eusBm)( 


(A^S) 


^’2  -  A;2,  Vs 
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k:^  -  kl 


The  vectors  ,  iV^  are  obtained  from  MmK  on  replacing  :^)  in 

the  expressions  for  An  etc.  in  (A-4)  by  -^YmikA.y).  Relations  (A-7)  are  also 

valid  when  Ymt{0,(p)  in  T~£,  etc.,  is  replaced  by  i,?)  cind  V^) 

Lm\  etc.,  is  replaced  by  Fr^’*(A:A,  <p).  In  (A-7),  P”‘  is  the  Legendre  associated 
polynomial. 


By  definition 


=  P^(cos^)(cosm(,5,sin77/.(,i>). 

The  associated  vector  harmonics  are  denoted  P^^  etc.  The  sanu'  a.pplie.s  to  cylin¬ 
drical  coordinates. 


The  recursion  relations  needed  for  our  computa 


tioiis  are: 


2£  + I 

2£-t-l(a^)  = - ze{x)  -  Zi_i{x), 

X 

2(.  +  l 

P^4.l(cos^)  =  — — —  cos6Pi{cos9)  — 


€  +  1 

9P£+i(cos^)  _  dPt^i{cos6) 

“  de 


•Pf_i(co.s6l), 


i  +  l 

—  {2C  -t-  l)(sin  0P(;{cos  0). 


(A -9) 


The  participating  scalar  eigenfunctions  are. 


xjo{x)  =  sin  I, 
jl{x)  =  sin  I  —  rcosr, 

X  72(3;)  =  (3  —  i"^)  sinx  —  3a:  cos  r, 
x^j^ix)  =  (15  —  6x^)  sin  x  +  (x^  —  15x)  cos  x, 
x^jiix.)  =  (105  —  45x^  +  x^)sinx  -f  (lOx^  —  105x)  cos  .r, 
x^js{x)  =  (945  —  420x^  +  15x'*)sinx  -|-  (— a;^  +  lOaa:'*  —  945a:)  cos  .r, 
a: A’6(a:)  =  (10,395  -  4725x2  +  210x'‘  -  X®)  sin X 
+  (— 21x^  4-  1260x^  —  10, 395x)  cos  x, 
x®j7(x)  =  (135,135  -  62,370x2  4-3150x4  -28x®) sin. r 

4-(x^  -  37Sx^  4-  17,325x2  -  135, 135.r)  cos  .r, 
x^jsix)  =  (2, 027, 025  -  945, 945x2  -h  51, 975x4  -  GOO.r"  +  )  sin  .r 

+  (36x^  -  6930x2  4-  270, 270x2  _  027^  025,;: )  cos  ,r. 
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(A -11) 


x^e^^h^^\x)  =  (-105x  +  lOx^)  +  ^■(105  -  45x2  + 
x®e‘^4^^(x)  =  (-x^  +  105x^  -  945x)  +  i(945  -  420x2  ^ 
x’’e‘®4^^(x)  =  (-21x^  +  1260x^  -  10, 395x) 

+  i(10, 395- 4725x2  +210x'^  -  x*'^), 
x®e‘^4^\x)  =  (x^  -  378x^  +  17, 325x2  -  135, 135x) 

+  i(135, 135  -  62,370x2  +  3150x'’  -  2Sx^), 
x^e^^hf\x)  =  (36x^  -6930x2+270,270x2  -  2,027,025x) 
+  i(2, 027, 025  -  945, 945x2  +  41, 975x'^ 

-  630x®  +  X®), 


Pq{cos6)  =  1, 

Pi  (cos  =  cos^, 

P2(cos^)  =  I  cos2  ^  ^  , 

5  3 

P3(cos  ^)  =  -  cos2  d  —  —  cos  6, 
nr  /i\  4  n  2/1  ^ 

P4(cos  0)  =  —  cos^  9  — —  cos"^  ^  +  o  ’ 

8  8  8 
*7(\  1  K 

P5(cos  ^)  =  —  cos®  6  — —  cos®  ^  "t"  "g"  ^5 

„  ,  231  6  z.  315  4  ^  105  2  z.  5 

P6(cos  6>)  =  —  cos®  ^  cos^  ^  ^  ~  16  ’ 

D  /  ziN  429  7  ^  693  5  ^  315  3  ^  35 

P7(cos  9)  =  — —  cos'  0  — —  cos®  6  H - —  cos®  9 - cos  9, 

16  16  16  16 

„  6435  sz,  12,012  g  z,  6930  4^ 

ft(«>s«)  =  — cos  9--^cos  «  +  ^cos  e 


Note  that 


1260  2  z.  35 

- cos®  9  H - 

128  128 


(A -12) 


(A -13) 


where  ai{x)  are  •polynomials  obeying  the  recursion  relation 

a;4.l(x)  =  {2i  +  l)a£(x)  -  x2af_](x), 

subjected  to  the  initial  conditions  ao(x)  =  0,  ai(x)  =  — x.  The  polynomials  t£(x) 


(A  -  14) 


obey  the  same  relation,  but  with  the  initial  conditions  bQ{x)  =  1,  bi(x)  =  1.  Also 

I 

ai{x)  +  ibiix)  =  J22-'^  (zx)^+i-^  (-4-15) 

Similarly,  if  we  put  in  (A-7) 

X  j£{x)  =  a^sinx  +  bicosx,  (.4  —  16) 

we  find  that  {a,  6}  obey  (A-14)  subjected  to  the  initial  conditions 

ao  =  1,  oi  =  1;  6o  =  0>  =  —X.  (A  —  17) 

In  the  far- field  we  may  assume  kpr  ^  1,  kgV  ^  1.  Consequently,  the  spherical 
Hankel  functions  appearing  in  (58)  may  be  approximated  l)y 

hf\x)  S  il^e-“+o(i), 


=  ^2  ^ 


e— +  0  4  , 


(A -18) 


Hence 


£o,o{kp^)  = 

’  g~zfcpr  ■ 

1  ^ 

L 

£o,2(kpr)  ^ 

'  ^-ikpr' 

[3  cos^  9 

kpr 

L  2 

^0,2i^sr)  ^ 

2  [  ksr 

sin  26eg . 

(A  -  19) 
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Appendix  B:  Quadrature  of  the  ‘Cavity  Integrals’ 


We  wish  to  evaluate  the  non-zero  spheroidal  cavity  coefficients  for  £  =  '2n, 
k  =  l  —  e^,ri  =  l  — 


{B  - 1) 


92n  =  ,  ;  X  /  f=^-P2n(cos0)sin  OcW.  {B-2) 

2n(2n  +1)7  —  n qqq2  q 

0  ^ 


Inserting  the  expansion 


P2n(c0S^)  =  C2n,2j  COS^-'  6,  (P  -  3) 

;=0 

where  C2n,2j  is  the  coefficient  of  the  Legendre  polynomial,  we  have 

7r/2 


/2n  =  (4n  +  1) 

i=0 
7r/2 


E^2„,2d/ 

n  ^  7 


cos 


2ie 


\/l  —  77  COs2  0 


sin  9(16 


COS^'^'^^  0 

fc  /  —  =  sin  /• , 


/ 

0 


\/l  —  ?7  cos^  6 


(B-4) 


_-(4t2-|-  Q  ^ 
72(2n-t-l)^^ 


n 

^ic2n,2j|  J  - 

i=l  ^  0 


COS 


2i0 


yjl  —  q  COs2  0 


sin 


7r/2 

-/: 


COS 


2i-h2  Q 


\/l  —  77  COs2  0 


sin 


} 


(P-5) 


It  is  thus  necessary  to  evaluate  only  a  single  integral-type  for  both  coefficients.  It 
is  readilv  shown  that 


7r/2 

/ 


cos 


2i  0 


\/l  —  77  cos2  0 


7r/2 

sin^d^  =  77“'^“^/^  J  cos^^  ^cl^,  [B  —  6) 


2.S 


/9  ■ 

where  cos(^0  =  and  cos^-^  is  given  in  (^).  Therefore, 


/2„  =  (477  -h  -  kq-^-'^^‘^l2j-^2m].  {B  -  7) 


_-(4n  +  l) 
92n  —  ,  1  \  ^ 


'■^^^/2i+2(?0)). 


:i?  -8) 


The  evaluation  of  the  corresponding  integrals  (37)-(3S)  for  the  circular  cylinder 
cavity  proceeds  along  the  same  lines.  Here  we  need,  in  addition,  the  exact  results 
(Magnus  et  al.,  1966) 


cos&P2n(cosff)sinffd0  =  — - — — — —  [-sin^  ^o^2n(cos0o  ) 

[Zn  —  i  ){Zn  +  Z) 

+  sin^  ^0  cos ^0-^271  (^°s  ^o)] , 


(B-9) 


sin^  0P2n(cos  0)  sin  ddff  =  —  sin^  OoP2n(cos  do ) 


(27?.  -  1)(27?  +  2) 


X  [sin^  ^o-P2n(cos  ^o)  —  cos  sin^  do )] , 


(B-10) 


(=— )• 

\  d  cos  d  J 

The  remaining  integrals  in  (37)-(38)  can  again  be  expressed  in  terms  of  l2j-  Al¬ 
together  we  find 


hn  =  (4"  +  i)[(5;7rTp;n:2jb”'«»^2„(coseo) 

n 

-sin'^  docosdoP2n(^os6o)]  -  '^C2n,2j{hj+2((^o)  “  hji^o)}  ^  -  H) 

i=0 

^2n  =  — j— ■  -  sin^  ^0-P2n(cos  Oq)  -  jr - iVo  T  doP2n(cOsdo) 

2n(2n  -b  1)  [  (2??  —  1)(27?  +  2)  ^ 


-  cos  dosin'^  do  P2n(cos  do)]  +'^2jc2n,2j{hj+2{So)  -  hjida)]  -{B  -12) 

i=i 
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Fig.  3 
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